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CO ■ In the primary visual cortex of primates and carnivores, functional architecture can be characterized 

I by maps of various stimulus features such as orientation preference (OP), ocular dominance (OD), and 

spatial frequency. It is a long-standing question in theoretical neuroscience whether the observed maps 
should be interpreted as optima of a speeifie energy functional that summarizes the design principles 
of cortical functional architecture. A rigorous evaluation of this optimization hypothesis is particularly 
demanded by recent evidence that the functional architecture of orientation columns precisely follows 
species invariant quantitative laws [T]. Because it would be desirable to infer the form of such an op- 
timization principle from the biological data, the optimization approach to explain cortical functional 
architecture raises the following questions: i) What are the genuine ground states of candidate energy 
funetionals and how can they be calculated with precision and rigor? ii) How do differences in candidate 
optimization principles impact on the predicted map structure and conversely what can be learned about 
an hypothetical underlying optimization principle from observations on map structure? iii) Is there a 
way to analyze the coordinated organization of cortical maps predicted by optimization principles in 



2 



general? To answer these questions we developed a general dynamical systems approach to the combined 
optimization of visual cortical maps of OP and another scalar feature such as OD or spatial frequency 
preference. From basic symmetry assumptions we obtain a comprehensive phenomenological classification 
of possible inter-map coupling energies and examine representative examples. We show that each individ- 
ual coupling energy leads to a different class of OP solutions with different correlations among the maps 
such that inferences about the optimization principle from map layout appear viable. We systematically 
assess whether quantitative laws resembling experimental observations can result from the coordinated 
optimization of orientation columns with other feature maps. 

Author Summary 

Neurons in the visual cortex form spatial representations or maps of several stimulus features. How 
arc different spatial representations of visual information coordinated in the brain? In this paper, we 
study the hypothesis that the coordinated organization of several visual cortical maps can be explained 
by joint optimization. Previous attempts to explain the spatial layout of functional maps in the visual 
cortex proposed specific optimization principles ad hoc. Here, we systematically analyze how optimiza- 
tion principles in a general class of models impact on the spatial layout of visual cortical maps. For 
each considered optimization principle we identify the corresponding optima and analyze their spatial 
layout. This directly demonstrates that by studying map layout and geometric inter-map correlations 
one can substantially constrain the underlying optimization principle. In particular, we study whether 
such optimization principles can lead to spatially complex patterns and to geometric correlations among 
cortical maps as observed in imaging experiments. 

Introduction 

Neurons in the primary visual cortex are selective to a multidimensional set of visual stimulus features, 
including visual field position, contour orientation, ocular dominance, direction of motion, and spatial 
frequency [2l|3]. In many mammals, these response properties form spatially complex, two-dimensional 
patterns called visual cortical maps [4H24|. The functional advantage of a two dimensional mapping of 
stimulus selectivities is currently unknown |25H27j . What determines the precise spatial organization of 
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these maps? It is a plausible hypothesis that natural selection should shape visual cortical maps to build 
efficient representations of visual information improving the 'fitness' of the organism. Cortical maps are 
therefore often viewed as optima of some cost function. For instance, it has been proposed that cortical 
maps optimize the cortical wiring length [28(129] or represent an optimal compromise between stimulus 
coverage and map continuity [30H43| . If map structure was largely genetically determined map struc- 
ture might be optimized through genetic variation and Darwinian selection on an evolutionary timescale. 
Optimization may, however, also occur during the ontogenetic maturation of the individual organism for 
instance by the activity-dependent refinement of neuronal circuits. If such an activity-dependent refine- 
ment of cortical architecture realizes an optimization strategy its outcome should be interpreted as the 
convergence towards a ground state of a specific energy functional. This hypothesized optimized func- 
tional, however, remains currently unknown. As several different functional maps coexist in the visual 
cortex candidate energy functionals are expected to reflect the multiple response properties of neurons in 
the visual cortex. In fact, consistent with the idea of joint optimization of different feature maps cortical 
maps are not independent of each other [5]|10lll8ll22ll41ll4iH47j . Various studies proposed a coordinated op- 
timization of different feature maps [5IIl[5^l551l551[571l551 - HTlH51H8H5n] . Coordinated optimization appears 
consistent with the observed distinct spatial relationships between different maps such as the tendency of 
iso-orientation lines to intersect OD borders perpendicularly or the preferential positioning of orientation 
pinwheels at locations of maximal eye dominance [HllinillSlllllllIlllllllS]- Specifically these geometric 
correlations have thus been proposed to indicate the optimization of a cost function given by a compro- 
mise between stimulus coverage and continuity [5^IMl[571l391l41lH5] , a conclusion that was questioned by 
Carreira-Perpinan and Goodhill [5T] . 

Visual cortical maps are often spatially complex patterns that contain defect structures such as point sin- 
gularities (pinwheels) [51 [TTl[5^[551l 531154] or line discontinuities (fractures) [T^llSn] and that never exactly 
repeat [T]1^ H2^1I5S] . It is conceivable that this spatial complexity arises from geometric frustration due to a 
coordinated optimization of multiple feature maps in which not all inter-map interactions can be simulta- 
neously satisfied |50[l57f[5n] . In many optimization models, however, the resulting map layout is spatially 
not complex or lacks some of the basic features such as topological defects [IHlEnilSZlinilini] ■ In other 
studies coordinated optimization was reported to preserve defects that would otherwise decay |50[l57j . 
An attempt to rigorously study the hypothesis that the structure of cortical maps is explained by an 
optimization process thus raises a number of questions: i) What are the genuine ground states of candi- 
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date energy functionals and how can they be calculated with precision and rigor? ii) How do differences 
in candidate optimization principles impact on the predicted map structure and conversely what can be 
learned about an hypothetical underlying optimization principle from observations on map structure? iii) 
Is there a way to analyze the coordinated organization of cortical maps predicted by optimization princi- 
ples in general? If theoretical neuroscience was able to answer these questions with greater confidence, the 
interpretation and explanation of visual cortical architecture could build on a more solid foundation than 
currently available. To start laying such a foundation, we examined how symmetry principles in general 
constrain the form of optimization models and developed a formalism for analyzing map optimization 
independent of the specific energy functional assumed. 

Minima of a given energy functional can be found by gradient descent which is naturally represented by a 
dynamical system describing a formal time evolution of the maps. Response properties in visual cortical 
maps are arranged in repetitive modules of a typical spatial length called hypcrcolumn. Optimization 
models that reproduce this typical length scale are therefore effectively pattern forming systems with a 
so-called 'cellular' or finite wavelength instability, see [55H65| . In the theory of pattern formation, it is 
well understood that symmetries play a crucial role |63H65j . Some symmetries are widely considered bi- 
ologically plausible for cortical maps, for instance the invariance under spatial translations and rotations 
or a global shift of orientation preference [SHIIS^IBSHTO] . In this paper we argue that such symmetries and 
an approach that utilizes the analogy between map optimization and pattern forming systems can open 
up a novel and systematic approach to the coordinated optimization of visual cortical representations. 
A recent study found strong evidence for a common design in the functional architecture of orientation 
columns [T|. Three species, galagos, ferrets, and tree shrews, widely separated in evolution of modern 
mammals, share an apparently universal set of quantitative properties. The average pinwheel density as 
well as the spatial organization of pinwheels within orientation hypercolumns, expressed in the statistics 
of nearest neighbors as well as the local variability of the pinwheel densities in cortical subregions ranging 
from 1 to 30 hypercolumns, are found to be virtually identical in the analyzed species. However, these 
quantities are different from random maps. Intriguingly, the average pinwheel density was found to be 
statistical indistinguishable from the mathematical constant tt up to a precision of 2%. Such apparently 
universal laws can be reproduced in relatively simple self-organization models if long-range neuronal in- 
teractions are dominant [TI ISWTT] . As pointed out by Kaschube and coworkers, these findings pose strong 
constraints on models of cortical functional architecture [T]. Many models exhibiting pinwheel annihi- 
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lation [Sni[S7] or pinwhecl crystallization [STJ[S21[Z2] were found to violate the experimentally observed 
layout rules. In [J it was shown that the common design is correctly predicted in models that were based 
only on intrinsic OP properties. Alternatively, however, it is conceivable that they result from geometric 
frustration due to inter-map interactions and joint optimization. In the current study we therefore in 
particular examined whether the coordinated optimization of the OP map and another feature map can 
reproduce the quantitative laws defining the common design. 

The presentation of our results is organized as follows. First we introduce a formalism to model the 
coordinated optimization of complex and real valued scalar fields. Complex valued fields can represent 
for instance orientation preference (OP) or direction preference maps [T3l[23]. Real valued fields may 
represent for instance ocular dominance (OD) [2], spatial frequency maps [T9l|44] or ON-OFF segrega- 
tion [73] • We construct several optimization models such that an independent optimization of each map 
in isolation results in a regular OP stripe pattern and, depending on the relative representations of the 
two eyes, OD patterns with a regular hexagonal or stripe layout. A model-free, symmetry-based analysis 
of potential optimization principles that couple the real and complex valued fields provides a compre- 
hensive classification and parametrization of conceivable coordinated optimization models and identifies 
representative forms of coupling energies. For analytical treatment of the optimization problem we adapt 
a perturbation method from pattern formation theory called weakly nonlinear analysis |63fl65l[71H77) . 
This method is applicable to models in which the spatial pattern of columns branches off continuously 
from an unselective homogeneous state. It reduces the dimensionality of the system and leads to am- 
plitude equations as an approximate description of the system near the symmetry breaking transition 
at which the homogeneous state becomes unstable. We identify a limit in which inter-map interactions 
that are formally always bidirectional become effectively unidirectional. In this limit, one can neglect the 
backreaction of the complex map on the layout of the co-evolving scalar feature map. We show how to 
treat low and higher order versions of inter-map coupling energies which enter at different order in the 
perturbative expansion. 

Second we apply the derived formalism by calculating optima of two representative low order examples of 
coordinated optimization models and examine how they impact on the resulting map layout. Two higher 
order optimization models are analyzed in Text SI. For concreteness and motivated by recent topical 
interest P T751I75] , we illustrate the coordinated optimization of visual cortical maps for the widely studied 
example of a complex OP map and a real feature map such as the OD map. OP maps are characterized by 
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pinwheels, regions in which columns preferring aU possible orientations are organized around a common 
center in a radial fashion [5^15^1501151] . In particular, we address the problem of pinwheel stability in OP 
maps [SOIED] and calculate the pinwheel densities predicted by different models. As shown previously, 
many theoretical models of visual cortical development and optimization fail to predict OP maps pos- 
sessing stable pinwheels |28ll50l[57ll61| . We show that in case of the low order energies, a strong inter- map 
coupling will typically lead to OP map suppression, causing the orientation selectivity of all neurons to 
vanish. For all considered optimization models, we identify stationary solutions of the resulting dynamics 
and mathematically demonstrate their stability. We further calculate phase diagrams as a function of the 
inter-map coupling strength and the amount of overrepresentation of certain stimuli of the co-evolving 
scalar feature map. We show that the optimization of any of the analyzed coupling energies can lead 
to spatially relatively complex patterns. Moreover, in case of OP maps, these patterns are typically 
pinwheel-rich. The phase diagrams, however, differ for each considered coupling energy, in particular 
leading to coupling energy specific ground states. We therefore thoroughly analyze the spatial layout of 
energetic ground states and in particular their geometric inter-map relationships. We find that none of 
the examined models reproduces the experimentally observed pinwheel density and spatially aperiodic 
arrangements. Our analysis identifies a seemingly general condition for interaction induced pinwheel-rich 
OP optima namely a substantial bias in the response properties of the co-evolving scalar feature map. 

Results 

1.1 Modeling the coordinated optimization of multiple maps 

We model the response properties of neuronal populations in the visual cortex by two-dimensional scalar 
order parameter fields which are either complex valued or real valued [5^115^ . A complex valued field z (x) 
can for instance describe OP or direction preference of a neuron located at position x. A real valued field 
o(x) can describe for instance OD or the spatial frequency preference. Although we introduce a model 
for the coordinated optimization of general real and complex valued order parameter fields we consider 
z(x) as the field of OP and o(x) as the field of OD throughout this article. In this case, the pattern of 
preferred stimulus orientation ^ is obtained by 

z9(x) = larg(z). (1) 
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The modulus |z(x)| is a measure of the selectivity at cortical location x. 

OP maps are characterized by so-called pinwheels, regions in which columns preferring all possible ori- 
entations arc organized around a common center in a radial fashion. The centers of pinwheels arc point 
discontinuities of the field t9(x) where the mean orientation preference of nearby columns changes by 90 
degrees. Pinwheels can be characterized by a topological charge q which indicates in particular whether 
the orientation preference increases clockwise or counterclockwise around the pinwheel center, 

where Ci is a closed curve around a single pinwheel center at x^. Since i9 is a cyclic variable in the interval 
[0,7r] and up to isolated points is a continuous function of x. qi can only have values 

where n is an integer number [83]. If its absolute value \qi\ = 1/2, each orientation is represented only 
once in the vicinity of a pinwheel center. In experiments, only pinwheels with a topological charge of 
±1/2 are observed, which arc simple zeros of the field 2:(x). 

OD maps can be described by a real valued two-dimensional field o(x), where o(x) < indicates ipsilateral 
eye dominance and o(x) > contralateral eye dominance of the neuron located at position x. The 
magnitude indicates the strength of the eye dominance and thus the zeros of the field correspond to the 
borders of OD. 

In this article, we view visual cortical maps as optima of some energy functional E. The time evolution 
of these maps can be described by the gradient descent of this energy functional. The field dynamics 
thus takes the form 



dtz{x,t) = F[z(x,t),o(x,t)] 

dto{x,t) = G[z(x,i),o(x,i)], (4) 

where F[z,o] and G[2,o] are nonlinear operators given by F[z,o] = — G[z,o] = The system 

then relaxes towards the minima of the energy E. The convergence of this dynamics towards an attractor 
is assumed to represent the process of maturation and optimization of the cortical circuitry. Various 
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biologically detailed models have been east to this form [MIISOIIM] . 

All visual cortieal maps are arranged in repetitive patterns of a typical wavelength A. We splitted the 
energy functional E into a part that ensures the emergence of such a typical wavelength for each map and 
into a part which describes the coupling among different maps. A well studied model reproducing the 
emergence of a typical wavelength by a pattern forming instability is of the Swift-Hohenberg type [64ll85] . 
Many other pattern forming systems occurring in different physical, chemical, and biological contexts (see 
for instance [74H77] ) have been cast into a dynamics of this type. Its dynamics in case of the OP map is 
of the form 

dtz{^,t)=Lzi^,t)-\z\^z, (5) 
with the linear Swift-Hohenberg operator 

L = r-{kl + A)\ (6) 

kc = 27r/A, and N[z{x,t)] a nonlinear operator. The energy functional of this dynamics is given by 

E=-J £x (^z(x)Lz(x) - ^k(x)l^) . (7) 

In Fourier representation, L is diagonal with the spectrum 

\{k) ^ r ^ {kl ^ k'f . (8) 

The spectrum exhibits a maximum at fc = fcc-For r < 0, all modes are damped since \{k) < 0, Vfc and 
only the homogeneous state 2:(x) = is stable. This is no longer the case for r > when modes on 
the critical circle k = kc acquire a positive growth rate and now start to grow, resulting in patterns 
with a typical wavelength A = l-n jk^. Thus, this model exhibits a supercritical bifurcation where the 
homogeneous state looses its stability and spatial modulations start to grow. 
The coupled dynamics we considered is of the form 

9tz(x, i) = z(x, i) - |zpz - — 

oz 

SU 

9to(x,<) = Loo(x,i)-o3- — +7, (9) 

do 
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where i{o,z} = ''{o,^} ^ [J^c {o z} ^" ' ^^'^ 7 is a constant. To account for the species differences in 
the wavelengths of the pattern we chose two typical wavelengths ~ 2Tr /k^.z and — 27r/fcc.o- The 
dynamics of z(x, t) and o(x, t) is coupled by interaction terms which can be derived from a coupling 
energy U. In the uncoupled case this dynamics leads to pinwheel free OP stripe patterns. 

1.2 Symmetries constrain inter-map coupling energies 

How many inter-map coupling energies U exist? Using a phenomenological approach the inclusion and 
exclusion of various terms has to be strictly justified. We did this by symmetry considerations. The 
constant 7 breaks the inversion symmetry o(x) = — o(x) of inputs from the ipsilateral (o(x) < 0) or 
contralateral (o(x) > 0) eye. Such an inversion symmetry breaking could also arise from quadratic terms 
such as o(x)^. In the methods section we detail how a constant shift in the field o(x) can eliminate the 
constant term and generate such a quadratic term. Including either a shift or a quadratic term thus 
already represents the most general case. The inter-map coupling energy U was assumed to be invariant 
under this inversion. Otherwise orientation selective neurons would, for an equal representation of the 
two eyes, develop different layouts to inputs from the left or the right eye. The primary visual cortex 
shows no anatomical indication that there are any prominent regions or directions parallel to the cortical 
layers |66] . Besides invariance under translations Tyz{x) = z(x — y) and rotations R^z(x) = z{R'^^x) of 
both maps we required that the dynamics should be invariant under orientation shifts 2;(x) — > e*''z(x). 
Note, that the assumption of shift symmetry is an idealization that uncouples the OP map from the 
map of visual space. BresslofF and coworkers have presented arguments that Euclidean symmetry that 
couples spatial locations to orientation shift represents a more plausible symmetry for visual cortical 
dynamics |67l[86|, see also [87|. The existence of orientation shift symmetry, however, is not an all or 
none question. Recent evidence in fact indicates that shift symmetry is only weakly broken in the spatial 
organization of orientation maps |88[l89j . A general coupling energy term can be expressed by integral 
operators which can be written as a Volterra series 

00 „ Mo u„+u-/2 u 

^ /]^d^Xio(xi) Y\ d^^i ^i^j) II (fxkz{-Xk)Ku{xi,...,y.u), (10) 

U=Ua+U^'^ i=l j=Ua + l k = U„+U^/2 + l 

with an w-th. order integral kernel A'„. Inversion symmetry and orientation shift symmetry require Uq to 
be even and that the number of fields z equals the number of fields z. The lowest order term, mediating 
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an interaction between the fields o and z is given by w = 4, = 2 i.e. 

y (i^xid^.T2d^a;3(i^a;4 o(xi)o(x2)2;(x3)z(x4) Jsr4(xi,X2,X3,X4) . (11) 

Next, we rewrite Eq. (|lip as an integral over an energy density \J . We use the invariance under translations 
to introduce new coordinates 





= (1/4) 




3 


yi 


— Xi '^ra 




— X2 — y^rn 




X3 X^7j^ . 



(12) 



This leads to 



£^4 = y y (fyiSy2(P'y?. o(yi + x„j)o(y2 + x,,„)z(y3 + x„) 
3 

z(x„i - ^(yj - x,„))A'(yi,y2,y3) 

i 

d^XmUiix^). (13) 



The kernel K may contain local and non-local contributions. Map interactions were assumed to be local. 
For local interactions the integral kernel is independent of the locations y^. We expanded both fields in 
a Taylor series around x,„ 

^(x,„ + yi) = z{xm) + Vz(x,„)yi + . . . , o(x,„ + y,;) = o(x,„) + Vo(x„)y,; + . . . (14) 

For a local energy density we could truncate this expansion at the first order in the derivatives. The 
energy density can thus be written 

U4{xm) = J d?yid^y2d^ys (o(x,„) + Vo(x,„)yi) (o(x™) + Vo(x™)y2) (15) 
(z(x,„) + Vz(x„)y3) ^z(x„,) - Vz(x,„) ^(y, - x,„)^ -?^(yi, y2, ys) • 
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Due to rotation symmetry this energy density should be invariant under a simultaneous rotation of both 
fields. From all possible eombinations of Eq. ([T5|) only those are invariant in whieh the gradients of the 
fields appear as scalar products. The energy density can thus be written as 

U4 = f{ci,C2,...,Cs) 

= f{o'^,z'^,zz,oz,Vo-Vo,Vz-Vz,Vz-Vz,Vo-Vz) , (16) 

where we suppress the argument x,„. All combinations Cj can also enter via their complex conjugate. 
The general expression for L/4 is therefore 

^ - E ^'1'^^^^ + E ^''^^'^ + E '^f ^^^^ ■ (17) 

From all possible combinations we selected those which are invariant under orientation shifts and eye 
inversions. This leads to 

Ui = ho'^ + l2\z\'^ + h{'^o^of +h\Vz^z\^ 

+h{Vz^zf + Iq{Vo-Vo)o^ + Z7(Vz-Vz)|z|^ 
+lsiVz-Vz)z'^ + Iq{Vz-Vz)z'^ 
+;io(Vo-Vz)oz + Zii(Vo-Vz)oz 
+/i2o2|zp + /i3|Vo-Vzp + luiVz-Vz)o^ 

+ii5(yo-yo)\z\^ + ii6{Vz-Vz)(yo-Vo) . (is) 

The energy densities with prefactor li to do not mediate a coupling between OD and OP fields and 
can be absorbed into the single field energy functionals. The densities with prefactors 1^ and Ig (also with 
ho and In) are complex and can occur only together with Ig ~ Ig {ho = hi) to be real. These energy 
densities, however, are not bounded from below as their real and imaginary parts can have arbitrary 
positive and negative values. The lowest order terms which are real and positive definite are thus given 
by 

Ui = /i2o2|zp + /i3|V-oVzp + Zi40^Vz-Vz + /i5Vo-Vo|zp + he (Vz-Vz) (Vo-Vo) . (19) 
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The next higher order energy terms arc given by 

Ue = o2|z|4 + |z|2o4 + o^Vz-Vz + . . . (20) 

Here the fields o(x) and z(x) enter with an unequal power. In the corresponding field equations these 
interaction terms enter either in the linear part or in the cubic nonlincarity. We will show in this article 
that interaction terms that enter in the linear part of the dynamics can lead to a suppression of the 
pattern and possibly to an instability of the pattern solution. Therefore we considered also higher order 
interaction terms. 

These higher order terms contain combinations of terms in Eq. (|19p and arc given by 

+ (Vz-Vzf (Vo-Vo)^ + o2|zp|Vo-Vzp + . . . (21) 

As we will show below examples of coupling energies 

U = ao2|z|2 + (3 |Vz-Vop + to'^\z\'^ + e |Vz-Vo|^ , (22) 

form a representative set that can be expected to reproduce experimentally observed map relationships. 
For this choice of energy the corresponding interaction terms in the dynamics Eq. ([9]) arc given by 

SU 

^ = Na[o,o, z] + Nglo^o, z] + NJo,o,o,o, z, z,^] + Nrlo, o,o,o, Z, Zj'z] 

oz 

= -ao^z + /3V{aVo) + e2V {\a\^aVo) ~2to^\z\'^z, 

^jj 

— — = iVQ[o, z, z] + iVfl[o, z, z] + iVc[o, o, o, z, z, z, z] + 7V,-[o, o, o, z, z, z, zl 

00 

= -ao|zp+/3V(aVz) + e2V(|a|2aVz) -2to^|z|'^ + c.c. (23) 

with a = Vz-Vo and c.c. denoting the complex conjugate. In general, all coupling energies in Ui,UQ, 
and C/g can occur in the dynamics and we restrict to those energies which are expected to reproduce 
the observed geometric relationships between OP and OD maps. It is important to note that with this 
restriction we did not miss any essential parts of the model. When using weakly nonlinear analysis the 
general form of the near threshold dynamics is insensitive to the used type of coupling energy and we 
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Figure 1. Pinwheel annihilation, preservation, and generation in numerical simulations for 
different strengths of inter-map coupling and OD bias, Tq = 0.2, = 0.02. Color code of OP map with 
zero contours of OD map superimposed. A j = 0, f3 — B j ~ 0, (3 — 0.3 C and D 7 = 0.15, [3 = 0.3. 
Initial conditions identical in A - C, T/ = lO"', kc,o = kc,z- 



therefore expect similar results also for the remaining coupling energies. 

Numerical simulations of the dynamics Eq. (jQ]), see [62l|9^, with the coupling energy Eq. ([22]) and 
a = e = r = arc shown in Fig. [TJ The initial conditions and final states arc shown for different 
bias terms 7 and intcr-map coupling strengths /?. We observed that for a substantial contralateral bias 
and above a critical intcr-map coupling pinwheels are preserved from random initial conditions or are 
generated if the initial condition is pinwheel free. Without a contralateral bias the final states were 
pinwheel free stripe solutions irrespective of the strength of the inter-map coupling. 



1.3 Calculating ground states by coupled amplitude equations 

We studied Eq. ^ with the low order inter-map coupling energies in Eq. (j22p using weakly nonlinear 
analysis. We therefore rewrite Eq. ^ as 



(9to(x,t) = roo(x,i) - i^o(x,t) + A^2,«[o,o] - A^3,„[o,o, o] - iV3^c[o, z,z] , 



(24) 
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where we shifted both hncar operators as — + L'^, Lq = I'o + L^. The constant term 7 in Eq. © 
is replaced by a quadratic interaction term N2.u[o,o\ = 70^ with 7 oc see Methods. The uncoupled 
nonlincarities arc given by iV3.„ = \z\'^z, N^^u = while iVs.c and A^a.c are the nonlincarities of the low 
order inter-map coupling energy Eq. (P^ . We study Eq. (|24p close to the pattern forming bifurcation 
where and are small. We therefore expand both control parameters in powers of the small expansion 
parameter /i 

rz = t^Tzi + /i^r^2 + ^J.^rz3 + ... 

To = M'^oi + ^Ji^ro2 + M^r-os + . . . . (25) 

Close to the bifurcation the fileds are small and thus nonlincarities are weak. We therefore expand both 
fields as 

0(x, t) = ^Ol (x, t) + 11^02 (x, t) + yU^Os (x, t) + . . . 

2:(x, t) = nzi (x, t) + n^Z2ix, t) + M^Z3(X, t) + ... (26) 
We further introduced a common slow timescale T = r^t and insert the expansions in Eq. p4| and get 
= 

(^-L^Z2 + TzlZl - Tzidrz^ 

+fi^ [~rz2dTZi + rz2Z\ + r^iZ2 - rz\dTZ2 - ^"za - A'"3,u[2i, ^i, zi]) 
(-iV3,c[^;i,oi,oi]) 



(27) 
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and 

= ^L"oi 

(-i"o2 + Toioi - r^idroi + Vm^oI + ^^^^^02 + ■ • •^2.m[oi, oi]^ 
+fi^ (^-r^2dT0i +ro20i + roi02 - r ^^18x02 - -£"03 - ^3,m[oi, oi, oi]^ 

'■■ (28) 

Wc consider amplitude equations up to third order as this is the order where the nonlinearity of the low 
order inter- map coupling energy enters first. For Eq. (|114p and Eq. ()115|) to be fulfilled each individual 
order in fi has to be zero. At linear order in /i we get the two homogeneous equations 

L°2i=0, L% = 0. (29) 

Thus zi and oi are elements of the kernel of and L'^. Both kernels contain linear combinations of 
modes with a wavevector on the corresponding critical circle i.e. 

n 

zi{x,T) = ^ (Aj-(r)e*'=^-^ + ylj-(r)e-'^^-^) 
3 

n 

oi(x,r) = ^ (Bj(T)e^^^^ + B,(T)e'**^^-^) , (30) 
j 

with the complex amplitudes Aj = ^je"^J, Bj = Bje^'^^ and kj = fcc,z (cos(j7r/n), sin(j7r/n)), kj = 
kc,o (cos(j7r/n), sin(j7r/n)). In view of the hexagonal or stripe layout of the OD pattern shown in Fig. [TJ 
7T, = 3 is an appropriate choice. Since in cat visual cortex the typical wavelength for OD and OP maps 
are approximately the same |561l91j i.e. fcc,o = ^c.z = the Fourier components of the emerging pattern 
are located on a common critical circle. To account for species differences we also analyzed models with 
detuned OP and OD wavelengths in part (II) of this study. 
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At second order in we get 

LPz2 + r^izi - r^idrzi = 

^"02 + roioi - r^iSrOi = 0. (31) 

As zi and oi are elements of the kernel rzi = Toi = 0. At third order, when applying the solvability 
condition (see Methods), we get 

rz2dTZi = rz2Zi- PcN3,u[zi,zi,zi]~ PcN3,c[zi,oi,oi] (32) 

rz2dTOi = ro20i - y/r^ PcN2,u[0l,0i] - Pc^3,u[oi, Oi, Oi] - Pc^3,c[0l,2;i,Zi] . 

We insert the leading order fields Eq. (|117p and obtain the amplitude equations 

rz2dTA, = rz2Ai - ^ h,j\Bj \^A, - hBfA,- - ^ 1 A^f - ^ f.^AjAj-A,- - 
3 j 3 

rz2dTB, = ro2B, - ^ /ly | - 2^B,+iB,+2 - I^J l'^* " ■ • ■ • (^3) 

3 3 

For simplicity we have written only the simplest intcr-map coupling terms. Depending on the configu- 
ration of active modes additional contributions may enter the amplitude equations. In addition, for the 
product-type coupling energy, there are coupling terms which contain the constant S, see Methods and 
Eq. pO)) . The coupling coefhcients are given by 
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From Eq. (|120[) we see that inter- map coupling has two effects on the modes of the OP pattern. First, 
inter-map coupling shifts the bifurcation point from to ^r^ — ■ This can cause a potential 

destabilization of pattern solutions for large inter-map coupling strength. Second, inter-map coupling 
introduces additional resonant interactions that for instance couple the modes Ai and their opposite 
modes Ai- . In case of A <C -B <C 1 the inter-map coupling terms in dynamics of the modes B are small. 
In this limit the dynamics of the modes B decouples from the modes A and we can use the uncoupled 
OD dynamics, see Methods. When we scale back to the fast time variable and set rz2 = r^, ro2 = ro we 
obtain 

dtA, = r,A,~Yl h^J\B,\^A, ~ hBjA,- - ^ g,, \A,f A, - ^ A,- A,- 
i 'J 3 

r.dtB, = roB,-Y,h^J\A,\^B,-2^oB^+{B,+2-Y.a^]\Bo?B,. (35) 

j 3 

The amplitude equations are truncated at third order. If pattern formation takes place far from threshold 
fifth order, seventh order, or even higher order corrections have to be considered. Moreover, at seventh 
order the nonlinearitics of the higher order inter-map coupling terms enter the expansion. A derivation 
of amplitude equation with higher order inter-map coupling energies is presented in Text SI. 

1.4 Interpretation of coupling energies 

Using symmetry considerations we derived inter-map coupling energies up to eighth order in the fields, 
see Eq. ([TO)) . Eq. ipO)) . and Eg. ([?!]) . Which of these various optimization principles could reproduce 
realistic inter-map relationships such as a uniform coverage of all stimulus features? We identified two 
types of optimization principles that can be expected to reproduce realistic inter-map relationships and 
good stimulus coverage. First, product-type coupling energies of the form U = o^"|zp", n — 1,2,.... 
These energies favor configurations in which regions of high gradients avoid each other and thus leading 
to high coverage. Second, gradient-type coupling energies of the form U = |Vo-Vzp", n = 1,2,.... 
In experimentally obtained maps, iso-orientation lines show the tendency to intersect the OD borders 
perpendicularly. Perpendicular intersection angles lead to high coverage as large changes of the field z in 
one direction lead to small changes of the field o in that direction. To see that the gradient-type coupling 
energy favors perpendicular intersection angles we decompose the complex field z{x.) into the selectivity 
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\z\ and the preferred orientation d. Wc obtain 



U = |Vz-Vo| 



2n 



(|Vo-V In |z|p + 4|Vo-Vi?p) 



(36) 



If the orientation selectivity is locaUy homogeneous, i.e. Vlnj^j « 0, then the energy is minimized if 
the direction of the iso-orientation lines (Vi?) is perpendicular to the OD borders. In our symmetry- 
based analysis we further identified terms that are expected to lead to the opposite behavior for instance 
mixture terms such asU = o'^'S/ z^'z. 

Pinwheels are prominent features in OP maps. We therefore also analyze how different optimization 
principles impact on the pinwheel positions with respect to the co-evolving feature maps. The product- 
type coupling energies are expected to favor pinwheels at OD extrema. Pinwheels are zeros of z and are 
thus expected to reduce this energy term. They will reduce energy the most when |o| is maximal which 
should repel pinwheels from OD borders, where o(x) is zero. Also the gradient-type coupling energy is 
expected to couple the OD pattern with the position of pinwheels. To see this we decompose the field z 
into its real and imaginary part 



At pinwheel centers the zero contours of Rez and Imz cross. Since there VRcz and Vlmz arc almost 
constant and not parallel the energy can be minimized only if |Vo| is small at the pinwheel centers, i.e. 
the extrema or saddle-points of o(x). 

From the previous considerations we assume all coupling coefficients of the energies to be positive. A 
negative coupling coefficient can be saturated by higher order inter-map coupling terms. In the following, 
we discuss the impact of the low order inter-map coupling energies on the resulting optima of the system 
using the derived amplitude equations. The corresponding analysis for higher order inter-map coupling 
energies is provided in Text SI. 

2.1 Optima of particular optimization principles: Low order coupling terms 

As indicated by numerical simulations and weakly nonlinear analysis of the uncoupled OD dynamics, see 
Methods, we discussed the influence of the OD stripe, hexagon, and constant solutions on the OP map 




(37) 
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using the coupled amplitude equations derived in the previous section. A potential baekreaction onto 
the dynamics of the OD map can be neglected if the modes Aj of the OP map are much smaller than 
the modes Bj of the OD map. This can be achieved if ^ r^. We first give a brief description of the 
uncoupled OP solutions. Next, we study the impact of the low order coupling energies in Eq. (|22p on 
these solutions. We demonstrate that these energies can lead to a complete suppression of orientation 
selectivity. In the uncoupled case there are for > two stable stationary solutions to the amplitude 
equations Eq. (??), namely OP stripes 

z,t = Ae<'^'-^+^), A^V^., (38) 

and OP rhombic solutions 

2 

z,h = -4 ^ (e''=^ ■^+'^^- + e"'^^- •^+'^.- ) , (39) 
i=i 

with + = 00, (f>2 + 4>2- = 00 + 71', 00 an arbitrary phase, and A = -^/r^/S ~ 0.447y^. In 
the uncoupled case the angle a = arccosfci • k2/k'^ between the Fourier modes is arbitrary. The stripe 
solutions are pinwheel free while the pinwheel density for the rhombic solutions varies as p = 4 sin a and 
thus < p < 4. For the rhombic solutions pinwheels are located on a regular lattice. We therefore refer 
to these and other pinwheel rich solutions which are spatially periodic as pinwheel crystals (PWC). In 
particular, we refer to pinwheel crystals with as rhombic spatial layout as rPWC solutions and pinwheel 
crystals with a hexagonal layout as hPWC solutions. Without inter-map coupling, the potential of the 
two solutions reads Vst = —\r1 < V^h = ~f ^z, thus the stripe solutions are always energetically preferred 
compared to rhombic solutions. 

In the following we study three scenarios in which inter-map coupling can lead to pinwheel stabilization. 
First, a deformation of the OP stripe solution can lead to the creation of pinwheels in this solution. 
Second, inter-map coupling can energetically prefer the (deformed) OP rhombic solutions compared to 
the stripe solutions. Finally, inter-map coupling can lead to the stabilization of new PWC solutions. 
For the low order interaction terms the amplitude equations are given by dtAi = —SV/SAi, dtBi = 
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—5V/5Bi with the potential 

j 

+2a5 {Ai'A^-B^i + A1A3-B2 + A1-A2B3 + A1-A3B2) 

+ E ffg^ 1^^ P + E a^^^^^^B, + E gS^A-A^- B-B, 

+ E glfA^j-Wj + E 9^?^^A- B^B, , (40) 

with the uncoupled contributions 

3 3 ,3 

3 3 

Vb - -r^El^.f + jE5.,|S.nS.f ■ (41) 
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The couphng coefficients read g-j^ = 2a + 2/3 cos'^{aij), g-J^ = 2q + /3 (l + cos^(Q;.y )), = 2a + 
/? (l + cos^(aij)), ^Ij^'' = a + /3, where is the angle between the wavevector fc,; and kj. 

2.2 Product-type energy ?7 = ao^|zp 

We first studied the impact of the low order product-type coupling energy. Here, the constant S{'y) enters 
explicitly in the amplitude equations, see Eq. (HUl) and Eq. (1551) . 

2.2.1 Stationary solutions and their stability 

In the case of OD stripes, see Methods, with Bi = B,B2 = B3 = we get the following amplitude 
equations 

dtAi {r, ~a6^ -2a\B\^)Ai-aB^Ai- + net. 

dtA2 = {r^.-ad'^ -2a\B\^) A2~2a5A^-~B + TLCt. 

dt A3 = ir, - a5^ - 2a\B\^) A3 - 2aM2- B + net. (42) 
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where net. refers to non inter-map eoupling terms — jylj — Y^^^^i fij^j^j- ^i- j resulting from 

the potential Va, see Eq. (|126p . The equations for the modes A^- are given by interehanging the modes 
Ai and A^- as well as interehanging the modes Bi and Bi. The OP stripe solution in ease of inter- map 
coupling is given by 

z = ^le^'^i-^+'^i) + ^ (43) 

with Ai =p^/^ /{2s/2B^a), Ax- = p/V^, and 

p ~ Tz — 2B^a — a5^ — y/ (r^ — aS'^){rz — a{AB'^ + S'^)) and the phase relation 

— (j>i- = 2ipi + TT. In the uncoupled ease (a = 0) they reduce to Ai- = and Ai = y/r\. With 
increasing inter-map coupling the amplitude Ai- grows and the solutions are transformed, reducing the 
representation of all but two preferred orientations. The parameter dependence of this solution is shown 
in Fig.[2]A for different values of the bias 7. Both amplitudes become identical at a = Vz/iAB"^ +(5^) with 



Ax=Ax-=\- . (44) 



This pattern solution finally vanishes at 

a, = rz/{B'' +8'') = 'drzlr,. (45) 

This existence border is thus independent of the OD bias 7. Above this eoupling strength only the trivial 
solution Aj = 0, Vj is stable. 

In addition to the OP stripe patterns there exist rhombic OP solutions, see Fig. [2^. These rhombic 
solutions are pinwheel rich with a pinwhcel density of p = 4sin7r/3 ~ 3.46 but are energetically not 
preferred compared to the stripe solutions, see Fig. The rhombic solutions in the uncoupled case 
Ai ~ Ai- — A2 — A2- , A3 — A^- = are transformed by inter- map coupling. The phase relations are 
given by 

(j)l+<f)X- =00, 01 - 01- = 2-01 + TT 

02 + 02- = 00 +71", 03 + 03- = 00 -I- TT 
Ipl + (f>2 + (f>3 = 00, (46) 
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where (j)o is an arbitrary phase. Stationary ampUtudes are given by 





a{B^ - (52) 
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^2 = A2- = 

with (/ = v/-360062ck2j2 ^ (_g^^ ^ 3662a + 9a(52)2. With increasing inter-map couphng strength a the 
amphtudes A2 = A2- are suppressed, see Fig. ([2])B. In addition, for nonzero bias 7, there is an increase 
of the amphtudes A3 = A3- . The amphtudes A2 and A3 collapse at a/r, = 3/(12^2 + 2QB5 + 352)_ a 
further increase of the inter-map coupling strength leads to a suppression of these amplitudes and finally 
to the OP stripe pattern where ^2 = A3 0. 

In the case the OD map is a constant, Eq. (j9ip . the amplitude equations simplify to 

dt A, = (r, - a5^) A, - ^ .9»j|A,f A, - ^ A,A,- A,- . (48) 

3 3 

Thus inter-map coupling in this case only renormalizes the bifurcation parameter and the energetic ground 
state is thus a stripe pattern with an inter-map coupling dependent reduction of the amplitudes 

Ai = v/rr~c^,A2 = A3 = 0. (49) 

Therefore at Uc = r^jS^ the stripe pattern ceases to exist and the only stable solution is the trivial one 
i.e. Ai = 0. In addition, there is the rhombic solution with the stationary amplitudes 



,A3 = 0. (50) 

In the case of OD hexagons Bi = 6e''^% Eq. ((90|) . the amplitude equations read 

dt A, = (r, - 6a62 - ocb") A, - aS^ (a.-c^"^' + 2A(,+i)-e*('^'+'^'+^) + 2A(,+2)- e^(^'+'^'+=)) 
-2aB'^ (A,+ie*('^'-'^'+i) A,+2e'W''-'^'+^)) 

-2am (A(,;+i)-e-''^'+^ ^ A(,+2)- e-^'^'+i) + net. , (51) 
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where the indices are cyclic i.e. i + 3 = i. These amphtude equations have stripc-Hke solutions as well 
as solutions with a rhombic layout of the form Ai = Ai- = A2 = A2- , A3 = .43- . For both solutions 
the stationary phases depend on the inter-map coupling strength a. In this case stationary solutions of 
Eq. ()51|) are calculated numerically using a Newton method and initial conditions close to these solutions. 
In contrast to the case of OD stripes and OD constant solutions the amplitude equations ((5T|) have an 
additional type of PWC solution which have uniform amplitudes, Aj ~ Ae^'^' . The dynamics of their 
phases is given by 

dt = 2A^ ^ sin + ~ - 0^- ) 

-B^a ^ (2 sin (0i - - tpi + ipj) + 2 sin ((j>i - (j)^- -ipi-ipj)) 
—B^a sin {(pi — (j)^- — 2ipi) 

-2SaB (sin (0^ - (f>(i+i)- + V'i+2) + sin (0^ - 0(^+2)- + V^i+i)) • (52) 

When solving the amplitude equations numerically we observe that the phase relations vary with the 
intcr-map coupling strength for non-uniform solutions. But for the uniform solution the phase relations 
are independent of the inter-map coupling strength. The phases of the uniform solution are determined 
up to a free phase ip which results from the orientation shift symmetry z — > 2; e"^ of Eq. (jH]). We therefore 
choose (f)i = ipi- As an ansatz for the uniform solutions we use 

(f>j = V'j + U - 1)2tt/3 + A (5j- 2 

= -V, + U - l)27r/3 + A (5,- 1 + ^,,3) , (53) 

where Sij is the Kronecker delta and A a constant parameter. Note, that z(x) cannot become real since 
(pj —'Pj-- The equation for the uniform amplitudes is then given by 

dtA^ r^A - ^A^ ~ AaB^A - aS^A + ABa [B - 25) cos A , (54) 



while the phase dynamics reads 

dt(t)j =-Ba{B- 26) sin A. 



(55) 
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The stationarity condition is fulfilled for an arbitrary S only if A = or A = tt. The corresponding 
amplitudes are given by solving the stationarity condition for the real part and read 

/r,-a(3B^ + 2BS + S^) ^ / ^ a {5B^ ^ 2BS + S^) 
Aa=o ^ Y 9 ' -^^^'^ = V 9 • ^^^^ 

We calculate the stability properties of all solutions by linear stability analysis considering perturbations 
of the amplitudes Aj ^ .4 + aj, Aj- A + a.j- and of the phases 4>j — > + Lpj, (pj- — > (pj- + tpj- . 
This leads to a perturbation matrix M. Amplitude and phase perturbations in general do not decouple. 
We calculated the eigenvalues of the perturbation M matrix numerically and checked the results by 
direct numerical simulation of the amplitude equations. In case of the uniform solutions Eq. (|134[) the 
perturbation matrix M is explicitly stated in Text S3. 

The stability of the A = Oand A = tt uniform solutions depends on the coupling strength a and on the 
sign of (B — 26). As the solution of B{-j) = 2(5(7) is given by 7 = 7* in the stability range of OD hexagons 
there is only one possible stable uniform solution, the A ~ n uniform solution. This solution ceases to 
exist at < a (bB^ — 2B6 + (5^). This existence border is in fact independent of the bias 7 and given by 

ac = 3r^/ro. (57) 

Thus the limit — > makes the uniform solution unstable for smaller and smaller coupling strengths. 
2.2.2 Bifurcation diagram 

The case of OP solutions when interacting with OD stripes is shown in Fig. [2]A,B. In case of OP stripes 
inter-map coupling suppresses the amplitude Ai of the stripe pattern while increasing the amplitude of the 
opposite mode Ai- . This transformation reduces the representation of all but two preferred orientations. 
When both amplitudes collapse the resulting OP map is selective only to two orthogonal orientations 
namely "d = (pi and d = (pi +7r/2. We refer to these unrealistic solutions as orientation scotoma solutions. 
The phase relations ensure that OD borders that run parallel to the OP stripes are located at the OP 
maxima and minima i.e. in the center of the orientation scotoma stripes. With increasing inter-map 
coupling, this orientation scotoma pattern is suppressed until finally all amplitudes are zero and only 
the homogeneous solution is stable. In case of OP rhombs inter-map coupling makes the rhombic pat- 
tern more stripe-like by reducing the amplitude A2 = A2- ■ The mode A^ — A^- which is zero in the 
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uncoupled case increases and finally collapses with the mode ^2- Increasing inter- map coupling more 
suppresses all but the two modes Ai = Ai-, leading again to the orientation scotoma stripe pattern. 
The parameter dependence of OP solutions when interacting with OD hexagons is shown in Fig. [2]C,D. 
OP stripe solutions became above a critical inter-map coupling strength unstable against PWC solutions. 
This critical coupling strength strongly depended on the OD bias. OP rhombic solutions also became 
unstable against PWC but for a lower coupling strength than the OP stripes. Thus there is at interme- 
diate coupling strength a bistability between stripe-like solutions and PWC solutions. The potential of 
the OP stripe and OP rhombic solutions is shown in Fig.[2f ,F. Stripes are energetically preferred in the 
uncoupled case as well as for small inter-map coupling strength for which they are stable. 
To summarize, stripe solutions were deformed but no pinwheels were created for this solution. The 
rhombic solutions were energetically not preferred for low inter-map coupling whereas for intermediate 
inter-map coupling these solutions lose pinwheels and became stripe solutions. Instead, additional pin- 
wheel rich solutions with a crystal layout became stable for intermediate inter-map coupling. For large 
inter-map coupling orientation selectivity was completely suppressed. 
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Figure 2. Stationary amplitudes with coupling energy U = a |zpo^, = 0.2. Solid (dashed) 
lines: stable (unstable) solutions to Eq. (|42|) (OD stripes) and Eq. (|5T|) (OD hexagons). A,B OD 
stripes, 7 = (blue), 7 = 7* (green), 7 = I.47* (orange). C,D OD hexagons, 7 = I.47* (blue), 7 ~ 87* 
(red). A,C Transition from OP stripe solutions, B,D Transition from OP rhombic solutions. E Potential 
for OP stripes (red) and OP rhombs (blue) interacting with OD stripes, 7 = 0. F Potential for OP 
stripes and OP rhombs interacting with OD hexagons. Arrows indicate corresponding lines in the phase 
diagram. Fig. ([3]). 
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Figure 3. Phase diagram with the coupUng energy U = ao^jzp, ro = 0.2, ^ Tq. Vertical 
lines: stability range of OD stripes, hexagons, and constant solution. Magenta (orange) line: transition 
of stripes (rhombs) to the orientation scotoma solution. Blue line: stability border for the A = tt 
uniform solution (hPWC). Green line: stability line of stripe- like solutions. Red line: pattern solutions 
ceases to exist, see Eq. (P5|) and Eq. ([57]). Blue region: stability region of hPWC, gray region: No 
pattern solution exists. 

2.2.3 Phase diagram 

The phase diagram as a function of the OD bias 7 and the inter-map coupling strength a for this coupling 
energy is shown in Fig. [31 When rescaling the inter- map coupling strength as ajr^ the phase diagram is 
independent of the bifurcation parameter of the OP map ■ Thus for fixed Tq ^ r r the phase diagram 
depends only on two control parameters 7/7* and a/vz- The phase diagram contains the stability borders 
of the uncoupled OD solutions 7*, 72 : 73 7 71 • They correspond to vertical lines, as they are independent of 
the inter-map coupling in the limit <C ^q. At 7 = 7* hexagons become stable. Stripe solutions become 
unstable at 7 = 73. At 7 73 the homogeneous solution becomes stable while at 7 = 74 hexagons lose 
their stability. In the units 7/7* the borders 72,73,74 vary slightly with Tq , sec Fig. |9l and are drawn 
here for Tq = 0.2. Colored lines correspond to the stability and existence borders of OP solutions. In 
the region of stable OD stripes the OP stripes run parallel to the OD stripes. With increasing inter-map 
coupling strength the orientation preference of all but two orthogonal orientations is suppressed. In the 
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Figure 4. Ipsi-center pinwheel crystal. A OP map, superimposed are tlie OD borders (gray), 90% 
ipsilateral eye dominance (black), and 90% contralateral eye dominance (white), = 0.2,7 = 87*. 
Dashed lines mark the miit cell of the regular pattern. B Distribution of preferred orientations. 

region of stable OD hexagons stripe-like OP solutions dominate for low inter-map coupling strength. 
Above a critical bias dependent coupling strength the A = tt uniform solution becomes stable (blue line). 
There is a region of bistability between stripe-like and uniform solutions until the stripe-like solutions lose 
their stability (orange line) . OP rhombic solutions lose their stability when the uniform solution becomes 
stable. Thus there is no bistability between OP rhombs and OP uniform solutions. As in the case of OD 
stripes the uniform solution becomes unstable at a = rz/{3B^). Also in the case of OD hexagons the 
inter-map coupling leads to a transition towards the trivial solution where there is no OP pattern at all. 
In case of the OD constant solution the OP map is a stripe solution. Pinwheel rich solutions thus occur 
only in the region of stable OD hexagons. In the following we discuss the properties of these solutions. 

2.2.4 Interaction induced pinwheel crystals 



The uniform solution Eq. (|134p with A = tt is illustrated in Fig. |4l For all stationary solutions the 
positions of the pinwheels are fixed by the OD map and there are no translational degrees of freedom. 
The unit cell (dashed line) contains 6 pinwheels which leads to a pinwheel density of p = 6 cos7r/6 ~ 5.2. 
Two of them are located at OD maxima (contra center) while one is located at an OD minimum (ipsi 
center). The remaining three pinwheels are located at OD saddle-points. Therefore, all pinwheels are 
located where the gradient of the OD map is zero. The pinwheel in the center of the OP hexagon is at the 
ipsilateral OD peak. Because these pinwheels organize most of the map while the others essentially only 
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match one OP hexagon to its neighbors we refer to this pinwheel crystal as the Ipsi- center pinwheel crystal. 
The iso-orientation Hnes intersect the OD borders (gray) exactly with a right angle. The intersection 
angles are, within the stability range of OD hexagons, independent of the bias 7. The remarkable property 
of perfect intersection angles cannot be deduced directly from the coupling energy term. The solution 
is symmetric under a combined rotation by 60° and an orientation shift by —60°. The symmetry of the 
pattern is reflected by the distribution of preferred orientations, see Fig. |4p. Although the pattern is 
selective to all orientations the six orientations + n^, n = 0, 5 are slightly overrepresented. 
To summarize, the low order product-type inter-map coupling leads in case of OD hexagons to a transition 
from pinwheel free stripe solutions towards pinwheel crystals. The design of the PWC is an example of 
an orientation hypercolumn dominated by one pinwheel. With increasing inter-map coupling the PWC 
solution is suppressed until only the homogeneous solution is stable. In case of OD stripes or the constant 
solution the OP solutions are pinwheel free stripe pattern. 

2.3 Gradient-type energy U = /3|Vz Vop 

When using a gradient-type inter-map coupling energy the interaction terms are independent of the OD 
shift 5. In this case, the coupling strength can be rescaled as (3B'^ and is therefore independent of the bias 
7. The bias in this case only determines the stability of OD stripes, hexagons or the constant solution. 

2.3.1 Stationary solutions and their stability 

A coupling to OD stripes is easy to analyze in the case of a gradient-type inter-map coupling. The 
energetically preferred solutions are OP stripes with the direction perpendicular to the OD stripes for 
which U — 0. This configuration corresponds to the Hubel and Wiesel Ice-cube model [2]. In addition 
there are rPWC solutions with the stationary amplitudes Ai ~ Ai- ~ \J {r^ + 2S^/3)/5, A2 = A2- = 
\J {rz — B'^P)/h, A^ = A3- = 0; and the stationary phases as in Eq. pS)) . Increasing inter-map coupling 
strength 13 leads to an increase of the amplitudes Ai = Ai- while decreasing the amplitudes A2 ~ A2- 
thus making the rhombic solution more stripe-like. 

In the case the OD map is a constant, Eq. (|9ip. the gradient- type inter- map coupling leaves the OP 
dynamics unaffected. The stationary states are therefore OP stripes with an arbitrary direction and 
rPWC solutions as in the uncoupled case. 
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In the case of OD hexagons the amphtudc equations read 




) 



+ net. . 



(58) 



Using Ai = Aic'^'^^ we obtain the phase equations 



Aidt (j>i 



^ AjAj-A^- sin + (/),- - <f)j - (j>j- ) 



-B'^P ^ f -Aj sin - (j>j -'(/'» + V'j) + -^Aj- sin - (j>j~ - tpi - ipj) 



-B^PA- sin((/., -4- -2i',) . 



(59) 



These amphtude equations have stripe-hke and rhombic solutions with inter-map couphng dependent 
phase relations. We therefore calculate their stationary phases and amplitudes numerically using a 
Newton method and initial conditions close to these solutions. Besides stripe-like and rhombic solutions 
these amplitude equations have uniform solutions. Again we find that the ansatz Eq. (|134p can satisfy 
the stationarity condition. The phase dynamics in this case reads 



As in the case of the product-type inter-map coupling energy stationary solutions are A = and A = tt 
with the stationary amplitudes 



We studied the stability properties of both stationary solutions by linear stability analysis where ampli- 
tude and phase perturbations in general do not decouple. The stability matrix of the uniform solutions 
is given in Text S3. The eigenvalues are calculated numerically. It turned out that the A = tt solution is 
unstable for /3 > while the A = solution becomes stable for (3 « O.OSr^/S^. The A = solution loses 
its stability above 




(60) 




(61) 



l3,B^^2r,/3. 



(62) 
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Figure 5. Stationary amplitudes with U = /?|Vz-Vop and OD hexagons. Solid (dashed) hnes: 
Stable (unstable) solutions to Eq. ([58|) . Transition from OP stripes towards the uniform solution (red), 
transition from OP rhombs towards the uniform solution (blue). Arrows indicate corresponding lines in 
the phase diagram, Fig. 

From thereon only the homogeneous solution Aj = is stable. 

2.3.2 Bifurcation diagram 

The course of the stationary amplitudes when interacting with OD hexagons is shown in Fig. [H The 
OP rhombic solution is almost unchanged by inter-map coupling but above a critical coupling strength 
the rhombs decay into a stripe-like solution. The amplitude of the OP stripe solution is suppressed by 
inter-map coupling and finally becomes unstable against the A = uniform solution. Thus for large 
inter-map coupling only the uniform solution is stable. A further increase in the inter-map coupling 
suppresses the amplitude of this uniform solution until finally only the homogeneous solution is stable. 

2.3.3 Phase diagram 

The phase diagram of this coupling energy is shown in Fig. [51 We rescaled the inter-map coupling 
strength as f3B'^/rz, where B is the stationary amplitude of the OD hexagons. The stability borders are 
then independent of the OD bias in the OD solutions and further independent of the bifurcation parameter 
rz- This simplifies the analysis since the OP solutions and their stability depend on 7 only indirect via 
the amplitudes B. In case of OD stripes or OD constant solution there is no pinwheel crystallization. 
Instead the OP solutions arc pinwheel free stripes. In case of OD hexagons hPWCs become stable above 
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Figure 6. Phase diagram with the coupUng energy U = /3|Vz-Vo|^, r, ^ r,,. Vertical black 
lines: stability range of OD stripes, hexagons, and constant solution. Blue line: stability border for the 
A = uniform solution. Green line: rhombic solutions become unstable. Orange line: stripe-like 
solutions become unstable. Red line: pattern solutions cease to exist, see Eq. (|6ip . Gray region: No 
pattern solution exists. 

/3 « 0.05rz/B^ (blue line). Rhombic OP patterns become unstable at (3B^/rz ~ 0.17 and decay into a 
stripe-like solution (green line). At /3B'^/rz w 0.36 these stripe-like solutions become unstable (orange 
line). Thus above ^S'^/r^ « 0.36 the hPWC is the only stable solution. At PB^/r^ = 2/3 the pattern 
solution ceases to exist (red line). 

2.3.4 Interaction induced pinwheel crystals 

The uniform solution Eq. (|134p . A = is illustrated in Fig. [T] This PWC contains only three pinwheels 
per unit cell leading to a pinwheel density of p = 3cos7r/6 ~ 2.6. Two of the three pinwheels arc located 
at maxima of the OD map (contra peak) while the remaining pinwheel is located at the minimum (ipsi 
peak) of the OD map. A remarkable property of this solution is that the pinwheel located at the OD 
minimum, carries a topological charge of 1 such that each orientation is represented twice around this 
pinwheel. Pinwheels of this kind have not yet been observed in experimentally recorded OP maps. This 
kind of vmiform solution corresponds to the structural pinwheel model by Braitenberg |92| . We therefore 
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Figure 7. The Braitenberg pinwheel crystal, A = uniform solution of Eq. (|134|) . A OP map, 
superimposed are the OD borders (gray), 90% ipsilateral eye dominance (black), and 90% contralateral 
eye dominance (white), ro = 0.2,7 = 37*. Dashed lines mark the unit cell of the regular pattern. B 
Distribution of preferred orientations. 

refer to this solution as the Braitenberg pinwheel crystal. 

The iso-orientation lines are again perfectly perpendicular to OD borders and this is independent of the 
bias 7. The solution is symmetric under a combined rotation by 120° and an orientation shift by — 27r/3. 
Further it is symmetric under a rotation by 180°. The pattern is selective to all orientations but the 
distribution of represented orientations is not uniform. The three orientations + n^, n = 0, ...,2 are 
overrepresented, see Fig. [7p. 

Overall this OP map is dominated by uniform regions around hyperbolic points. In contrast to the ipsi 
center PWC all pinwheels in this OP map organize a roughly similar fraction of the cortical surface. 

Discussion 

Summary of results 

In this study we presented a symmetry-based analysis of models formalizing that visual cortical archi- 
tecture is shaped by the coordinated optimization of different functional maps. In particular, we focused 
on the question of whether and how different optimization principles specifically impact on the spatial 
layout of functional columns in the primary visual cortex. We identified different representative candidate 
optimization principles. We developed a dynamical systems approach for analyzing the simultaneous op- 
timization of interacting maps and examined how their layout is infiuenced by coordinated optimization. 
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In particular, wc found that intcr-map coupling can stabilize pinwheel-rich layouts even if pinwheels are 
intrinsically unstable in the weak coupling limit. We calculated and analyzed the stability properties of 
solutions forming spatially regular layouts with pinwheels arranged in a crystalline array. We analyzed 
the structure of these pinwheel crystals in terms of their stability properties, spatial layout, and geomet- 
ric inter-map relationships. For all models, we calculated phase diagrams showing the stability of the 
pinwheel crystals depending on the OD bias and the inter-map coupling strength. Although differing in 
detail and exhibiting distinct pinwheel crystal phases for strong coupling, the phase diagrams exhibited 
many commonalities in their structure. These include the general fact that the hexagonal PWC phase is 
preceded by a phase of rhombic PWCs and that the range of OD biases over which pinwheel crystalliza- 
tion occurs is confined to the stability region of OD patch solutions. 
Comparison to previous work 

Our analytical calculations of attractor and ground states close a fundamental gap in the theory of visual 
cortical architecture and its development. They rigorously establish that models of interacting OP and 
OD maps in principle offer a solution to the problem of pinwheel stability [501170] . This problem and 
other aspects of the influence of OD segregation on OP maps have previously been studied in a series 
of models such as elastic net models [31l31linillIl[Sni[SH] , self-organizing map models [HTllMlllSlllHlllS] , 
spin-like Hamiltonian models [STUHS], spectral filter models [M], correlation based models [nSllHSl, and 
evolving field models |97| . Several of these simulation studies found a higher number of pinwheels per hy- 
percolumn if the OP map is influenced by strong OD segregation compared to the OP layout in isolation 
or the influence of weak OD segregation [SOII^IW] . In such models, large gradients of OP and OD avoid 
each other [351135]. As a result, pinwheel centers tend to be located at centers of OD columns as seen in 
experiments |18[|4411^51l47[l95] . By this mechanism, pinwheels are spatially trapped and pinwheel anni- 
hilation can be reduced [50] . Moreover, many models appear capable of reproducing realistic geometric 
inter-map relationships such as perpendicular intersection angles between OD borders and iso-orientation 
lines [HHlinilinS]. Tanaka et al. reported from numerical simulations that the relative positioning of ori- 
entation pinwheels and OD columns was dependent on model parameters [SS]. Informative as they were, 
almost all of these previous studies entirely relied on simulation methodologies that do not easily permit 
to assess the progress and convergence of solutions. Whether the reported patterns were attractors or 
just snapshots of transient states and whether the solutions would further develop towards pinwheel-free 
solutions or other states thus remained unclear. Moreover, in almost all previous models, a continuous 
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variation of the intcr-map coupling strength was not possible which makes it hard to disentangle the 
contribution of inter-map interactions from intrinsic mechanisms. The only prior simulation study of a 
coordinated optimization model that tracked the number of pinwheels during the optimization process 
did not provide evidence that pinwhecl annihilation could be stopped but only reported a modest reduc- 
tion in annihilation efficiency |50| . From this perspective, the prior evidence for coordination induced 
pinwhecl stabilization appears relatively limited. Our analytical results leave no room to doubt that map 
interactions can stabilize an intrinsically unstable pinwhecl dynamics. They also reveale that interaction 
of orientation preference with a stripe pattern of OD is per se not capable of stabilizing pinwheels. 
The mathematical structure of interaction models 

Independent of its predictions, our study clarifies the general mathematical structure of interaction dom- 
inated optimization models. To the best of our knowledge our study for the first time describes an 
analytical approach for examining the solutions of coordinated optimization models for OP and OD 
maps. Our symmetry-based phenomenological analysis of conceivable coupling terms provides a general 
classification and parametrization of biologically plausible coupling terms. To achieve this we mapped the 
optimization problem to a dynamical systems problem which allows for a perturbation expansion of fixed 
points, local minima, and optima. Using weakly nonlinear analysis, we derived amplitude equations as an 
approximate description near the symmetry breaking transition. We identified a limit in which inter-map 
coupling becomes effectively unidirectional enabling the use of the uncoupled OD patterns. We studied 
fixed points and calculated their stability properties for different types of inter-map coupling energies. 
This analysis revealed a fundamental difference between high and low order coupling energies. For the 
low order versions of these energies, a strong inter-map coupling typically leads to OP map suppression, 
causing the orientation selectivity of all neurons to vanish. In contrast, the higher order variants of the 
coupling energies do generally not cause map suppression but only infiuence pattern selection, see Text 
SI. We did not consider an interaction with the retinotopic map. Experimental results on geometric 
relationships between the retinotopic map and the OP map are ambiguous. In case of ferret visual cortex 
high gradient regions of both maps avoid each other f41| . In case of cat, however, high gradient re- 
gions overlap |17j . Such positive correlations cannot be easily treated with dimension reduction models, 
see |98| . It is noteworthy that our phenomenological analysis identified coupling terms that could induce 
an attraction of high gradient regions. Such terms contain the gradient of only one field and can thus be 
considered as a mixture of the gradient and the product-type energy. 
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Conditions for pinwheel stabilization 

Our results indicate that a patchy layout of a second visual map interacting with the OP map is im- 
portant for the effectiveness of pinwheel stabilization by inter-map coupling. Such a patchy layout can 
be easily induced by an asymmetry in the representation of the corresponding stimulus feature such as 
eye dominance or spatial frequency preference. In spatial frequency maps, for instance, low spatial fre- 
quency patches tend to form islands in a sea of high spatial frequency preference |44| . In our model, the 
patchy layout results from the overall dominance of one eye. In this case, OD domains form a system of 
hexagonal patches rather than stripes enabling the capture and stabilization of pinwheels by inter-map 
coupling. The results from all previous models did not support the view that OD stripes are capable 
of stabilizing pinwheels [50l|93l[97j . Our analysis shows that OD stripes are indeed not able to stabilize 
pinwheels, a result that appears to be independent of the specific type of map interaction. In line with 
this, several other theoretical studies, using numerical simulations [50119311971 . indicated that more banded 
OD patterns lead to less pinwheel rich OP maps. For instance, in simulations using an elastic net model, 
the average pinwheel density of OP maps interacting with a patchy OD layout was reported substantially 
higher (about 2.5 pinwheels per hypercolumn) than for OP maps interacting with a more stripe-like OD 
layout (about 2 pinwheels per hypercolumn) [50] . 

Experimental evidence for pinwheel stabilization by inter-map coupling 

Several lines of biological evidence appear to support the picture of interaction induced pinwheel sta- 
bilization. Supporting the notion that pinwheels might be stabilized by the interaction with patchy 
OD columns, visual cortex is indeed dominated by one eye in early postnatal development and has a 
pronounced patchy layout of OD domains [SSlllOlj . Further support for the potential relevance of this 
picture comes from experiments in which the OD map was artificially removed resulting apparently in a 
significantly smoother OP map |43j . In this context it is noteworthy that macaque visual cortex appears 
to exhibit all three fundamental solutions of our model for OD maps: stripes, hexagons, and a monocular 
solution, which are stable depending on the OD bias. In the visual cortex of macaque monkeys, all three 
types of patterns are found near the transition to the monocular segment, see [101| and Fig. ([5]). Here, 
OD domains form bands in the binocular region and a system of ipsilateral eye patches at the transition 
zone to the monocular region where the contralateral eye gradually becomes more dominant. If pinwheel 
stability depends on a geometric coupling to the system of OD columns one predicts systematic differences 
in pinwheel density between these three zones of macaque primary visual cortex. Because OD columns 
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in the binocular region of macaque visual cortex are predominantly arranged in systems of OD stripes 
our analysis also indicates that pinwheels in these regions are either stabilized by other patchy columnar 
systems or intrinsically stable. 

The geometry of interaction induced pinwheel crystals 

One important general observation from our results is that map organization was often not inferable 
by simple qualitative considerations on the energy functional. The organization of interaction induced 
hexagonal pinwheel crystals reveals that the relation between coupling energy and resulting map structure 
is quite complex and often counter intuitive. We analyzed the stationary patterns with respect to inter- 
section angles and pinwheel positions. In all models, intersection angles of iso-orientation lines and OD 
borders have a tendency towards perpendicular angles whether the energy term mathematically depends 
on this angle, as for the gradient-type energies, or not, as for the product-type energies. Intersection angle 
statistics thus are not a very sensitive indicator of the type of interaction optimized. Mathematically, 
these phenomena result from the complex interplay between the single map energies and the interaction 
energies. In case of the low order gradient-type inter-map coupling energy all pinwheels are located at 
OD extrema, as expected from the used coupling energy. For other analyzed coupling energies, however, 
the remaining pinwheels are located either at OD saddle-points (low order product-type energy) or near 
OD borders (higher order gradient-type energy) , in contrast to the expection that OD extrema should be 
energetically preferred. Remarkably, such correlations, which arc expected from the gradient-type cou- 
pling energies, occur also in the case of the product-type energies. Remarkably, in ease of product type 
energies pinwheels are located at OD saddle-points, which is not expected per se and presumably result 
from the periodic layout of OP and OD maps. Correlations between pinwheels and OD saddle-points 
have not yet been studied quantitatively in experiments and may thus provide valuable information on 
the principles shaping cortical functional architecture. 
How informative is map structure? 

Our results demonstrate that, although distinct types of coupling energies can leave distinguishing sig- 
natures in the structure of maps shaped by interaction (as the OP map in our example) , drawing precise 
conclusions about the coordinated optimization principle from observed map structures is not possible 
for the analyzed models. In the past numerous studies have attempted to identify signatures of coordi- 
nated optimization in the layout of visual cortical maps and to infer the validity of specific optimization 
models from aspects of their coordinated geometry [551[gSll571 l5gWT | Hg i li5H50 1 [57 1 [55 1 [g5] . It was, however. 
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never clarified theoretically in which respect and to which degree map layout and geometrical factors of 
inter-map relations are informative with respect to an underlying optimization principle. Because our 
analysis provides complete information of the detailed relation between map geometry and optimization 
principle for the different models our results enable to critically assess whether different choices of energy 
functionals specifically impact on the predicted map structure and conversely what can be learned about 
the underlying optimization principle from observations of map structures. 

We examined the impact of different interaction energies on the structure of local minima and ground 
states of models for the coordinated optimization of a complex and a real scalar feature map such as OP 
and OD maps. The models were constructed such that in the absence of interactions, the maps reorga- 
nized into simple stripe or blob pattern. In particular, the complex scalar map without interactions would 
form a periodic stripe pattern without any phase singularity. In all models, increasing the strength of 
interactions could eventually stabilize qualitatively different more complex and biologically more realistic 
patterns containing pinwheels that can become the energetic ground states for strong enough inter-map 
interactions. The way in which this happens provides fundamental insights into the relationships between 
map structure and energy functionals in optimization models for visual cortical functional architecture. 
Our results demonstrate that the structure of maps shaped by inter-map interactions is in principle in- 
formative about the type of coupling energy. The organization of the complex scalar map that optimizes 
the joined energy functional was in general different for all different types of coupling terms examined. 
We identified a class of hPWC solutions which become stable for large inter-map coupling. This class 
depends on a single parameter which is specific to the used inter-map coupling energy. Furthermore, as 
shown in Text SI, pinwheel positions in rPWCs, tracked while increasing inter-map coupling strength, 
were different for different coupling terms examined and thus could in principle serve as a trace of the 
underlying optimization principle. This demonstrates that, although pinwheel stabilization is not re- 
stricted to a particular choice of the interaction term, each analyzed phase diagram is specific to the used 
coupling energy. In particular, in the strong coupling regime substantial information can be obtained 
from a detailed inspection of solutions. 

In the case of the product-type coupling energies, the resulting phase diagrams are relatively complex as 
stationary solutions and stability borders depend on the magnitude of the OD bias. Here, even quan- 
titative values of model parameters can in principle be constrained by analysis of the map layout. In 
contrast, for the gradient-type coupling energies, the bias dependence can be absorbed into the coupling 
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strength and only selects the stationary OD pattern. This leads to relatively simple phase diagrams. For 
these models map layout is thus uninformative of quantitative model parameters. We identified several 
biologically very implausible OP patterns. In the case of the product-type energies, we found orientation 
scotoma solutions which are selective to only two preferred orientations. In the case of the low order 
gradient-type energy, we found OP patterns containing pinwheels with a topological charge of 1 which 
have not yet been observed in experiments. If the relevant terms in the coupling energy could be deter- 
mined by other means, the parameter regions in which these patterns occur could be used to constrain 
model parameters by theoretical bounds. 

The information provided by map structure overall appears qualitative rather than quantitative. In both 
low order inter-map coupling energies (and the gradient-type higher order coupling energy, see Text SI), 
hPWC patterns resulting from strong interactions were fixed, not exhibiting any substantial dependence 
on the precise choice of interaction coefficient. In principle, the spatial organization of stimulus prefer- 
ences in a map is an infinite dimensional object that could sensitively depend in distinct ways to a large 
number of model parameters. It is thus not a trivial property that this structure often gives essentially 
no information about the value of coupling constants in our models. The situation, however, is reversed 
when considering the structure of rPWCs. These solutions exist and are stable although energetically not 
favored in the absence of inter-map interactions. Some of their pinwheel positions continuously depend on 
the strength of inter-map interactions. These solutions and their parameter dependence nevertheless are 
also largely uninformative about the nature of the interaction energy. This results from the fact that rP- 
WCs are fundamentally uncoupled system solutions that are only modified by the inter-map interaction. 
As pointed out before, preferentially orthogonal intersection angles between iso-orientation lines and OD 
borders appear to be a general feature of coordinated optimization models in the strong coupling regime. 
Although the detailed form of the intersection angle histogram is solution and thus model specific, our 
analysis does not corroborate attempts to use this feature to support specific optimization principles, 
see also [5nil02iri03j . The stabilization of pinwheel crystals for strong inter-map coupling appears to be 
universal and provides per se no specific information about the underlying optimization principle. In fact, 
the general structure of the amplitude equations is universal and only the coupling coefficients change 
when changing the coupling energy. It is thus expected that also for other coupling energies, respecting 
the proposed set of symmetries, PWC solutions can become stable for large enough inter-map coupling. 
Conclusions 



40 



Our analysis conclusively demonstrates that OD segregation can stabilize pinwheels and induce pinwheel- 
rich optima in models for the coordinated optimization of OP and OD maps when pinwheels are intrin- 
sically unstable in the uncoupled dynamics of the OP map. This allows to systematically assess the 
possibility that inter-map coupling might be the mechanism of pinwheel stabilization in the visual cor- 
tex. The analytical approach developed here is independent of details of specific optimization principles 
and thus allowed to systematically analyze how different optimization principles impact on map layout. 
Moreover, our analysis clarifies to which extend the observation of the layout in physiological maps can 
provide information about optimization principles shaping visual cortical organization. 
The common design observed in experimental OP maps [1] is, however, not reproduced by the optima 
of the analyzed optimization principles. Whether this is a consequence of the applied weakly nonlinear 
analysis or of the low number of optimized feature maps or should be considered a generic feature of 
coordinated optimization models will be examined in part (II) of this study [90) . In part (II) we comple- 
ment our analytical studies by numerical simulations of the full field dynamics. Such simulations allow 
to study the rearrangement of maps during the optimization process, to study the timescales on which 
optimization is expected to take place, and to lift many of the mathematical assumptions employed by 
the above analysis. In particular, we concentrate on the higher order inter-map coupling energies for 
which the derived amplitude equations involved several simplifying conditions, see Text SI. 



We studied the intersection angles between iso-orientation lines and OD borders. The intersection angle 
of an OD border with an iso-orientation contour q;(x) is given by 



where x denotes the position of the OD zero-contour lines. A continuous expression for the OP gradient 
is given by Vt9 = ImVz/z. We calculated the frequency of intersection angles in the range [0,7r/2]. In 
this way those parts of the maps are emphasized from which the most significant information about the 
intersection angles can be obtained [18]. These are the regions where the OP gradient is high and thus 
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every intersection angle receives a statistical weight according to |Vi9|. For an alternative method see |10[ . 
The transition from OD stripes to OD blobs 

We studied how the emerging OD map depends on the overall eye dominance. To this end we mapped the 
uncoupled OD dynamics to a Swift-Hohenberg equation containing a quadratic interaction term instead 
of a constant bias. This allowed for the use of weakly nonlinear analysis to derive amplitude equations 
as an approximate description of the shifted OD dynamics near the bifurcation point. We identified the 
stationary solutions and studied their stability properties. Finally, we derived expressions for the fraction 
of contralateral eye dominance for the stable solutions. 

Mapping to a dynamics with a quadratic term 

Here we describe how to map the Swift-Hohenberg equation 

9to(x,t) = Lo(x,t) -o(x,t)^ +7, (64) 

to one with a quadratic interaction term. To eliminate the constant term wc shift the field by a constant 
amount o(x, t) = o(x, t) + S. This changes the linear and nonlinear terms as 

Lo — > Lo—{k^ — ro)5 

-d^ + 3Sd^ + 3d^d + . (65) 

We define the new parameters r,, = r,, — and 7 = —3(5. This leads to the new dynamics 

dtd^¥od- (fc2 + A)^5 + 7o2-o3-53- (fc^^^-r„)(5 + 7. (66) 

The condition that the constant part is zero is thus given by 



- (5^ - (fc^ - ro) J + 7 = 



(67) 
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For To < 1 the real solution to Eq. ([S7)) is given by 



^ 32VI' (68) 



1/3 



with /3 = 1^—277 + Y 108 (ro — fee) + 7297^J . For small 7 this formula is approximated as 

S-^TT^ 7'7Tr^^ + 37'7Tr^^ + --- (69) 

The uncoupled OD dynamics we consider in the following is therefore given by 

9(0 = foO-(A:2 + A)^o + 7o2-o3^ (70) 

This equation has been extensively studied in pattern formation literature [76] . 
Amplitude equations for OD patterns 

We studied Eq. (j70p using weakly nonlinear analysis. This method leads to amplitude equation as an 
approximate description of the full field dynamics Eq. (j70p near the bifurcation point Fq = 0. We 
summarize the derivation of the amplitude equations for the OD dynamics which is of the form 

dt o(x, t) = L o(x, t) + 7V2 [o, o] - N-i [o, o, o] , (71) 

with the linear operator L = — {k^ ^ + A)'^. In this section we use for simplicity the variables (o, Tq, 7) 
instead of {o,fo,"f)- The derivation is performed for general quadratic and cubic nonlinearities but are 
specified later according to Eq. ^ as A^a [0,0,0] = and iV2[o, o] =70^. For the calculations in the 
following, it is useful to separate Tq from the linear operator 

L^r, + L\ (72) 

therefore the largest eigenvalue of is zero. The amplitude of the field o(x, t) is assumed to be small 
near the onset = and thus the nonlinearities are small. We therefore expand both the field o(x, t) 
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and the control parameter in powers of a small expansion parameter fi as 

o(x, t)^fiOi (x, t) + n^02 (x, t) +n^03{x,t) + ... , (73) 

and 

To = f^ri + /i^rz + + . . . (74) 

The dynamics at the critical point Tq = becomes arbitrarily slow since the intrinsic timescale t = 
diverges at the critical point. To compensate we introduce a rescaled time scale T as 

T = rot, dt^ To Ot . (75) 

In order for all terms in Eq. (j7ip to be of the same order the quadratic interaction term N2 must be 
small. We therefore rescale N2 as ^/rZN2- This preserves the nature of the bifurcation compared to the 
case A^2 = 0. 

We insert the expansion Eq. ([73| and Eq. ([74| in the dynamics Eq. ([7T|) and get 

= ^i"oi 

+ (^-L^02 - ridroi + rioi + v^/iri + i^flr2 + . . . A^2[oi,oi]^ 
+ (^-Z^os + ri (02 - 0x02) + ^2 (oi - droi) - Nj,[oi,oi,oi]^ 

'■■ (76) 



We sort and collect all terms in order of their power in /i. The equation can be fulfilled for /i > only if 
each of these terms is zero. We therefore solve the equation order by order. In the leading order we get 
the homogeneous equation 

L°oi=0. (77) 

Thus oi is an element of the kernel of . The kernel contains linear combinations of modes with 
wavevector kj on the critical circle \kj\ — kc.o- At this level any of such wavevectors is possible. We 
choose 

n n 

«i ^^Bj{T)e'''^-^ + ^Bj{T)e-'''^-^ , (78) 
j 3 
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where the wavevectors arc chosen to be equally spaced kj = kc.o (cos(j7r/n), sin(j7r/n)) and the complex 
amplitudes Bj = Bje'^'^^ . The homogeneous equation leaves the amplitudes Bj undetermined. These 
amplitudes are fixed by the higher order equations. Besides the leading order homogeneous equation we 
get inhomogeneous equations of the form 

i°o,„ ^ Fm (79) 

To solve this inhomogeneous equation we first apply a solvability condition. We thus apply the Fredholm 
Alternative theorem to Eq. ([75]) . Since the operator LP is self-adjoint LP — LP\ the equation is solvable 
if and only if Fm is orthogonal to the kernel of LP i.e. 

(F„,5)=0, VL"6 = (80) 

The orthogonality to the kernel can be expressed by a projection operator Pc onto the kernel and the 
condition (F, o) = can be rewritten as PcF = 0. 
At second order we get 

L%2 = ri (oi - droi) . (81) 

Applying the solvability condition Eq. (|M| we see that this equation can be fulfilled only for ri =0. At 
third order we get 

L°03 = r2 (oi - dTOi)+N2[oi,oi] - N3[oi,oi,oi] . (82) 

The parameter r2 sets the scale in which oi is measured and we can set r2 — I. We apply the solvability 
condition and get 

dr oi = 01 + PcN2 [oi , oi] - A^a [oi , oi , oi] . (83) 
We insert our ansatz Eq. ([78| which leads to the amplitude equations at third order 

OtB, = B, + P, BjBke-'^^^N2[e'^'^, e'''"^] - P,J2 BjBkBie-^'^^^ N^le'^'^ , e'''"^ , e'^'^\ , (84) 

0,k j,k 

where Pi is the projection operator onto the subspace {e''^'^} of the kernel. Pi picks out all combinations 
of the modes which have their wavevector equal to ki . In our case the three active modes form a so called 
triad resonance ki + k2 + ^3 = 0. The quadratic coupling terms which are resonant to the mode Bi are 
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therefore given by 

B^B^e-'^^^ (N2[e-'^^^, e"*^^^] + iYsfe"*'^^^, e"*^^^]) . (85) 

Resonant contributions from the cubic nonhnearity resuh from terms of the form \Bj\'^Bi. Their couphng 
coefficients are given by 

and 

= iV3[e''^'^, e'*^-*, e"*'^-^] + iV3[e*'^-^, e"''^'^, e'^'""] + iV3[e-'^'^, e'^^^ , e'*^'*] . (87) 

When specifying the nonhnearities Eq. (|9]) the couphng coefficients are given by Ijij = 6,^,^ = 3. Finally, 
the amplitude equations (here in the shifted variables (^0,7) are given by 

dtBi = ¥oBi - 3|Si|2Bi - 6 (IB2P + IS3P) Si + 27B2B3 , (88) 

where we scaled back to the original time variable t. Equations for B2 and B3 are given by cyclic 
permutation of the indices. 

Stationary solutions 

The amplitude equations ((88|) have three types of stationary solutions, namely OD stripes 

Ost (x) = 2Bst cos (.T + V) + S, (89) 

with Bst = \/f/3, hexagons 

3 

o^,., (x) = Bh,, + ^-^^ + (90) 

with the resonance condition Yl^j^j — ^ and iJ/ie^r ~ — 7/I5 + \J (7/15)^ + r/15. Finally, there is a 
homogeneous solution with spatially constant eye dominance 

Oc(x) = 5. (91) 
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Figure 8. OD patterns. A Stationary amplitudes of the OD dynamics. The course oi Bst{j) Eq. ([55)1 
(blue), Bhexil) Eq. ([90]) (red), and of ^(7) Eq. (|68)) (green) for Tq = 0.2. The sohitions are plotted in 
solid lines within their stability ranges. B OD map of macaque monkey. Adapted from [101| . C Details 
of B with stripe-like, patchy, and homogeneous layout. 

The spatial average of all solutions is (o(x)) = S. The course of Bst, Bhex, and of 5(7) is shown in Fig. [8] 
Linear stability analysis for OD patterns 

We decomposed the amplitude equations (j88p into the real and imaginary parts. From the imaginary 
part we get the phase equation 

dtipi = -2jsm{ilji+ij2 + ^p3) , (92) 

and equations for ip2i4'3 by cyclic permutation of the indices. The stationary phases are given by tpi + 
"02 + V'3 = {0,7r}. The phase equation can be derived from the potential V[ip] = — 27Cos(V-'i + "02 + V'a)- 
We see that the solution tpi + ip2 + tpa = is stable for 7 > 0(7 < 0) and the solution ipi + '4'2 + tps — it 
is stable for 7 < 0(7 > 0). 

We calculate the stability borders of the OD stripe, hexagon, and constant solution in the uncoupled case. 
This treatment follows [76j . In case of stripes the three modes of the amplitude equations are perturbed 
as 

Bi^Bst + h, B2^b2, Bs^bs, (93) 
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assuming small perturbations fei, &2, and 63. This leads to the linear equations dfb ~ Mb with the stability 
matrix 



/ 



M 



7 -Wit 

r-&Bl 2jBst 

2^B,t j 

The corresponding eigenvalues are given by 



V 



(94) 



(95) 



This leads to the two borders for the stripe stability 



r = 0, = 37 ■ 



In terms of the original variables ^0,7 the borders are given by (0 < Tq < 1) 
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(3 - 2ro) . _ (15 - 14r„) 



33/2 



72 



153/2 



To derive the stability borders for the hexagon solution Ohex{'^) we perturb the amplitudes as 



(96) 



(97) 



Bi Bhex + ^1 , B2 ^ Bhex + ^2 , S3 — > Bhex + ^3 



(98) 



The stability matrix is then given by 
/ 

M = 



-21Bhex + r -l2Biex ~ "^Bhexl 

12^2^, - 2Bhexl -21Bhex + r 

y -12BL. - 2i?,e.7 -12i?2,.-2i?,.e.7 



-12i?L.-2B„.e.7 



-12i?L.-2B,e.7 
~2lBhex + r J 



(99) 



and the corresponding eigenvalues are given by 



A = (-455,2 +r- AByri, -9BI + r + 2Bhj, -9BI + r + 2Bhl) 



(100) 
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Figure 9. The uncoupled OD dynamics. A Phase diagram of the OD model Eq. (|64| . Dashed 
hnes: stabihty border of hexagon solutions, solid line: stability border of stripe solution, gray regions: 
stability region of constant solution B Percentage of neurons dominated by the contralateral eye plotted 
for the three stationary solutions. Circles: numerically obtained values, solid lines: Cst and Chex- 



The stability borders for the hexagon solution are given by 



1 ^2 ~ 16-2 

'15^' "^y^ ■ 



(101) 



In terms of the original variables we finally get 

(12 - 7ro) ^oVl 



24 



7 



(51 - 50r„) 
513/2 



(102) 



The phase diagram of this model is depicted in Fig.[51<\. It shows the stability borders 7*, 72 1 73: and 74 
for the three solutions obtained by linear stability analysis. Without a bias term the OD map is either 
constant, for ro < 0, or has a stripe layout, for > 0. For positive Tq and increasing bias term there are 
two transition regions, first a transition region from stripes to hexagons and second a transition region 
from hexagons to the constant solution. 

The spatial layout of the OD hexagons consists of hexagonal arrays of ipsilateral eye dominance blobs in 
a sea of contralateral eye dominance, see Fig. (SJA. 
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Contralateral eye fraction 

To compare the obtained solutions with physiological OD maps we quantified the fraction of neurons 
selective to the contralateral eye inputs. For stripe and hexagon solutions we thus calculated the fr-action 
of contralateral eye dominated territory Cgt and Chex- In case of stripes this is a purely one-dimensional 
problem. The zeros of the field are given by 

o,t {x) = 2Bst cos (x + 7/)) + 5 = , (103) 



with the solution 

T = a.rccos I 

As the field has a periodicity of tt the area fraction is given by 



X = arccos (77^ ) • (104) 



Cst = arccos (^^^^ /'^ ■ (105) 

In case of hexagons we observe that the territory of negative o(x) values is approximately a circular area. 
We obtain the fraction of negative o(x) values by relating this area to the area of the whole hexagonal 
lattice. In case of hexagons the field is given by 



Ohe 



;(x) = 2Bhex J2 {^1^ + ^j) + ^ ■ (106) 



As an approximation we project the field onto the a;-axis and choose for simplicity -ipj = 0, Vj. The field 
has its maximum at the origin Ohexi^, 0) ^ 6 + 6. The projection leads to 

f{x) = 2Bhex (cosx + 2cos(l/2)a;) +(5. (107) 

The zeros f{xi) = are located at 



.,=2arccos(^^-l + y3+^jj . (108) 

The circular area of positive Ohex (x) values is now given by Ac ~ Tixf . The periodicity of the hexagonal 
pattern is given by f{x2)s=o ~ min(/)5=o = ~3. This leads to X2 = 47r/3. The area of the hexagon is 
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therefore given by A^ex = 3x2^/5/2. The contra fraction is finally given by 

The course of the fractions Cst and Chex is shown in Fig. [9^. At the border 7 = 7*, where hexagons 
become stable Chex ~ 65.4%. At the border 7 = 74, where hexagons loose stability Chex ~ 95.2%. Both 
quantities are independent of Tq. Wc confirmed our results by direct numerical calculation of the fraction 
of positive Ohex{^) pixel values. Deviations from the result Eq. (|109p are small. For 7/7* ~ 1 the zeros of 
Eq. (jl07p are not that well approximated with a circular shape and the projection described above leads 
to the small deviations which decrease with increasing bias 7. 
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Text SI: Supporting Information for 

Coordinated optimization of visual cortical maps 
(I) Symmetry-based analysis 

3 Introduction 

In case of the low order intcr-map coupling energies strong inter-map coupling leads to a suppression of OP 
selectivity. This suppressive effect can be avoided by restricting coupling strengths. One aim of this article 
is to test different optimization principles and potentially rule out some optimization principles. When 
comparing our results from different optimization principles to biological data such parameter tuning 
reduces the practicability. In this supporting information we complement our study using the high order 
inter-map coupling energies. We show that in this case a suppression of OP selectivity cannot occur. 
We derive coupled amplitude equations which, however, involve several mathematical assumptions. A 
systematic treatment as it is shown in the main article would imply that low order and higher order inter- 
map coupling energies are in general non-zero. Low order energy terms would enter at third order in the 
expansion and higher order corrections could potentially alter the stability properties. In addition, higher 
order inter-map coupling energies can affect the stability of patterns. In the following, we assume all low 
order inter-map coupling energies to be zero and that contributions entering the amplitude equations at 
higher orders can be neglected. The obtained results are confirmed numerically in part (II) of this study. 

4 Coupled amplitude equations: Higher order terms 

We studied the coupled Swift-Hohenberg equations 



9to(x, t) = roo(x,t) - L° o(x,t) -I- iV2,„[o, o] - iV3,„[o, o, o] - 7V7,c[o, o, o, z, z, z, z] , (110) 
with the higher order inter-map coupling energies 



z(x,i) - iV3,„[z,z,z] - Nt,c[z,z,z, o,o, o,o] 



^7 = ro4|z|4 + e|Vz•Vo|^ 



(111) 
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using weakly nonlinear analysis. We study Eq. (|110p close to the pattern forming bifurcation where and 
To are small. We therefore expand both control parameters in powers of the small expansion parameter 
M 

Tz = l^Tzl + H^rz2 + ^I.^rz3 + ■■■ 

To = ^roi + Ai^ro2 + Ai^ro3 + . . . . (112) 

Close to the bifurcation the fileds are small and thus nonlinearities are weak. We therefore expand both 
fields as 

0(X, t) = ^Oi (X, t) + I? 02 (X, t) + /i^03 (x, t) + . . . 

z(x,t) ^izi{^,t) + ^? Z2{yi,t) + Z3,{-K,t) + ... (113) 
We further introduced a common slow timescale T — r^t and insert the expansions in Eq. (|110p and get 

= 

+^? (^-L^Z2 + Tzizi - rzidrzij 

+^-1^ (^-rz2dTZi + rz2Zi + r.iZ2 - rzidTZ2 - L^z^ - A^3,„[zi, zi, zi]^ 
(^-L°zj + rz2Z5 + rz4z3 + r^gZi H h N3^u[z5, Zi,zi]j 

(-iV7,c[2;i,Zl,Zl,Oi,Oi,Oi,Oi]) 



(114) 
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and 



= ^lL°Ol 

-rz2dTOi + ro20i + roi02 - rzidT02 - L^o^ - iV3,„[oi, oi, oi] 



(^-i/°07 + 7'o205 + fo403 + TogOi H ^3,« [05 , Ol , Oi] - 7V2,„ [01,05] " • ■ ■) 

+^J.'^ (|-^7x[oi , Ol , Ol , Zl, Zl , Zl , Zl]^ 



(115) 



We consider amplitude equations up to seventh order as this is the order where the nonUnearity of the 
higher order couphng energy enters first. For Eq. (|114p and Eq. pi5p to be fulfilled each individual order 
in ^ has to be zero. At linear order in fi we get the two homogeneous equations 

L°zi = 0, L%^0. (116) 

Thus zi and oi are elements of the kernel of ^^^d L^. Both kernels contain linear combinations of 
modes with a wavevector on the critical circle i.e. 



-J- )e 



zi(x,r) = Y,(^Af\T)e''''-^ + Af]{T)e 
j 

n 

oi{x,T) = 5^(Bf)(r)e'^~^-">;Bf^(T)e-^^^-^) , (117) 

j 

with the complex amplitudes Aj"'^' = Aje^'^\ Bj^'^ — Bje^'^^ and kj = fc^^^ (cos(j7r/n), sin(j7r/n)), k'^ = 
kc,o (cos(j'7r/n), sin(j7r/n)). In view of the hexagonal or stripe layout of the OD pattern shown in Fig. [l] 
71 = 3 is an appropriate choice. Since in cat visual cortex the typical wavelength for OD and OP maps 
are approximately the same |561l91j i.e. kc^o = ^c.z = the Fourier components of the emerging pattern 
are located on a common critical circle. To account for species differences we also analyzed models with 
detuned OP and OD wavelengths in part (II) of this study. 
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At second order in we get 

LPz2 + r^izi - r^idrzi = 

i"o2+roiOi-r,iaTOi = 0. (118) 

As zi and oi are elements of the kernel rzi = Toi = 0. At third order, when applying the solvability 
condition (see Methods), we get 

r^z^TOi = ro20i - ^7^Fc^2,«[oi,oi] - Pc^3,«[oi,Oi,oi] . (119) 
We insert the leading order fields Eq. (|117p and obtain the amplitude equations 

j 3 

r..drJ^' ^ -2./T^;^i,-Yj^^^^ ■ {120) 

j 

These uncoupled amplitude equations obtain corrections at higher order. There arc fifth order, seventh 
order and even higher order corrections to the uncoupled amplitude equations. In addition, at seventh 
order enters the nonlinearity of the higher order inter-map coupling energies. The amplitude equations 
up to seventh order arc thus derived from 

rz2dT zi = rz2Zi- PcN3^u[zi,zi,zi] 

rz2dT Z3 = rz2Z3 P7V3,„[zi, zi, zjj] (121) 

rz2dTZ5 = rz2Z5 PcA^s.uN, ^i, ^] - AA'"7,c[zi, Zl, ^l, Ol, Oi, Ol, Ol] , 

and corresponding equations for the fields oi, 03, and 05. The field zi is given in Eq. (jll7p and its ampli- 
tudes and B^i) arc determined at third order. The field Z3 contains contributions from modes off the 
critical circle Z3 o//, \koff\ 7^ and on the critical circle i.e. 23 = Z3 off+Yl ( A[^\T)e^''^'^ + A''^] {T)e^^''^'^ 
Its amplitude A^^^^ are determined at fifth order. The field Z5 also contains contributions from modes off 
the critical circle 25,0// and on the critical circle i.e. = ^5,0// + E (Af\T)e'''^-^ + A^^] {T)e-'''^-A . 
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Its amplitude A^^^ are determined at seventh order. This leads to a series of amplitude equations 



j j 

r..dr^' = r^.Af Y.^.^^^ (122) 

j jlk 

which are solved order by order. We set — Tz and = fo and reseale to the fast time. This leads to 
dt^' = rz^' Y.a^A^'^ (123) 

j jlk 

We can combine the amplitude equations up to seventh order by introducing the amplitudes Aj = 
A^P + Af' + A^J'^ and B, = s'lp + + . This leads to the amplitude equations 

dt Ai = TzAi ^ ^ 9ij\Aj I Ai ^ ^ fijAjAj- A^- 

-Y,K,iu\A,\^\Bi\^\Bu?A,-... 
Jlk 

dtB, = roB,-2B,+{B,+2-Y.9^^BJ?AB, 

3 

-Y,h^3lk\B,\Ml\Mk?B,- ... . (124) 

Jlk 

For simplicity we have written only the simplest inter-map coupling terms. Depending on the configu- 
ration of active modes additional contributions may enter the amplitude equations. In addition, for the 
product-type coupling energy, there are coupling terms which contain the constant 5, see Methods. In 
case of A ^ -B <C 1 the inter-map coupling terms in dynamics of the modes B are small. In this limit 
the dynamics of the modes B decouples from the modes A and we can use the uncoupled OD dynamics, 
see Methods. In the following, we use the effective inter-map coupling strength eB^ (and tB"^). 
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5 Higher order inter-map coupling energies 

5.1 Optima of particular optimization principles: Higher order coupling 



In this article we demonstrated that the low order coupling terms can lead to a complete suppression of 
OP selectivity i.e. vanishing magnitude of the order parameter \z\. As the coupling terms are effectively 
linear they not only influence pattern selection but also whether there is a pattern at all. This is in 
general not the case for higher order coupling energies using the amplitude equations Eq. (I124p . In this 
case the coupling is an effective cubic interaction term and complete selectivity suppression is impossible. 
Moreover, as in the low-order energy case, we could identify the limit <C Tq in which the backreaction 
onto the OD map formally becomes negligible. The potential is of the form 



terms 



V = Va + Vb + 





j,l,k u,v,w 



E ^-p A.AiAkMBuB, Sj:^ 




j,l,k u,v 



j,l,k u 



(125) 



where Sij denotes the Kronecker delta and the uncoupled contributions 



Va 





(126) 



Amplitude equations can be derived from the potential by dtAi = —SV/6Ai. We have not written 
terms involving the modes Aj- or Bj. The complete amplitude equations involving all modes and the 
corresponding coupling coefficients are given in Text S2. As for the low order coupling energies terms 
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involving the constant 5 depend only on the coupling coefficient of the product-type energy r. In the fol- 
lowing we specify the amplitude equations for negligible backreaction where B = Bjiex, B = Bgt, or B = 0. 

5.2 Product-type energy U = r o^\z\^ 

First, we studied the higher order product-type inter- map coupling energy TJ = t o^\z\^ . As for the lower 
order version of this coupling energy the shift (5(7) explicitly enters the amplitude equations resulting in 
a rather complex parameter dependence, see Eq. (68) in the Methods section. 

5.2.1 Stationary solutions and their stability 

In the case of OD stripes the amplitude equations of OP modes read 

^B^Al-M- J2 ^uA,Aa[{SS''B + 2ASB^B)Sj:^^j:^_j:^^^ (127) 

+ {SS'B + 2ASBB') + 8^i?^%„+£„_£„._2fcO 

dt A, = r,A, - ^ (.gg^ \A,\'A, + g!^^A,A^-A,- ) 
j 

-glf^A^A.-A.^g^'^A^AiA- - \l1gf A\--A^- - Xl^g^^^A^A^- 

+fl('°^,5 

"fc„+fc„-fe„,-/ci+A:2 /c„ + fe„-/c„,,2fei+A;2 ^ ^'U "fe„ +fe„ - A;„ ,-2fei +fe2 i ' 
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where 6i,j denotes the Kronecker delta and 



9lP= l + S^ + 12S^\B\^+6\B\\ pgi=2gW, 

slf = 9^^, = 12S^B^ + 8B^, 

0, 4J, = 2 + I2|i3|4 + 24<52|B|2 + 2J4^ 

9^= 0, 9^}^^ = 125'B^ + 8B^B^ 

5(5) = 12<52# + SBb\ g^:l = &\B\' + &5\B?, 9^'^ = 2gil\ 

5^1^ = + IBB^S, = 2541^ , 

g'il ^ 4i?j3 + xB^Bb, g'll^^ = 25!?,) , 

y«u — o , 9uv^u ~ ^9uu , 

5«« - bB b , 9uv^u - 29uu ■ 

The equation for the mode is given by interchanging the modes A2 and in Eq. I\\27\ . The equations 
for the modes A^- are given by interchanging the modes Ai and Aj- and interchanging the modes Bi 
and Bi. 

In this ease, at low inter-map coupling the OP stripes given by 

z = Aie^^'-'^'^''"') - Ai-e-<^'-^+'^^-) , (128) 

with (pi — (pi- = 2'ipi + TT run parallel to the OD stripes. Their stationary amplitudes are given by 

A\ = (^u — V — \/ — 2uv + v'^ — 16w^^ x/32'w 
Ai- = x/2, (129) 

with X = [u - V + yjv? - 2uv + - I6w^) /{uv - - 8w^), u = 2 + ViB^r + 2AB'^5'^t + 26^t, v = 
{QB'^ + 12B^6'^ + (5'^)r, w — {2B^ + 3S'^)tB'^. The parameter dependence of these stripe solutions is shown 
in Fig. [TOA. 

At large inter-map coupling the attractor states of the OP map consist of a stripe pattern containing 
only two preferred orientations, namely 1} = (pi and "i? = 0i + 7r/2. The zero contour lines of the OD map 
are along the maximum amplitude of orientation preference minimizing the energy term. 
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In addition there are rhombic solutions 

z = Aie'^'''-^+'^'''> + Ai-e-'^'''-^-'l''+'''> + A2e'^''^-^+'^'''> + y^2- e"'^^"'^"'''^^ , (130) 

which exist also in the uncoupled case, see Fig. HOB. However, these rhombic solutions are energetically 
not favored compared to stripe solutions, sec Fig. HOC . The inclusion of the intcr-map coupling makes 
these rhombic solution even more stripe-like. 
In case of a OD constant solution the amplitude equations read 

dtA, = r,A, - J2 9v\Aj\^A, - ^ f.^AjA^-A,- , (131) 
i j 

with gu = 1 + (S^T, gij =2 + 2S^t and /y = 2 + 2S'^t. Inter-map coupling thus leads to a renormalization 
of the uncoupled interaction terms. Stationary solutions are stripes with the amplitude 

and rhombic solutions with the stationary phases + — 4>2 — (t>2- = ^nd the stationary amplitudes 

Ai=Ai-=A2= A2- = yjr.jih + 5(54r) . (133) 

In the case of OD hexagons wc identify, in addition to stripe-like and rhombic solutions, uniform solutions 
Ai = A. When solving the amplitude equations numerically we have seen that the phase relations vary 
with the intcr-map coupling strength t for non-uniform solutions. But for the uniform solution the phase 
relations are independent of the inter-map coupling strength. We use the ansatz for uniform solutions 

v4, J ■ — ■ j — - — " ^ J - — " 1 ^ 2 ^ 3 
(l^j = + U - l)27r/3 + A (5j- 2 
<^j- = -V-, + (j - l)27r/3 + A (,5,- 1 + J,, 3) , (134) 
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where di_j is the Kroneckcr delta and A a constant parameter. This leads to the stationarity condition 



Four types of stationary solutions exist namely the A = 0, A = tt, which we already observed in case of 
the low order energies, and the solutions 



which depends on B and S and thus on the bias 7. The course of Eq. (|136p as a function of 7 is shown 
in Fig. [TTB. Stationary amplitudes for these solutions are given by 



3t (3/t + 33^4 + 56^3(5 + bOB'^S^ + 16BS^ + 3S^) 

!i (i^ij) 

r (9/r + 483i34 - 50463^ + 24662^2 _ 4gg^3 + 954) V y 

r,{13B^ - 8BS + 6S^) 

3t (41166t + 704B^St - 37664(5^ + 32B^6^t + 2S^{9 - 7(5V) - 39B^{3S^t - 1) + 86(5(5(54r - 3)) 



We study the stability properties of OP stripe-like, rhombic and uniform solutions using linear stability 
analysis. The eigenvalues of the stability matrix , see Text S2, are calculated numerically. Linear stability 
analysis shows that for r > the A ~ solution is unstable for all bias values. The stability region of 
the A — IT solution and the solution Eq. (|136p is bias dependent. The bias dependent solution Eq. (|136l) 
is stable for 7 > 7* and 7 < 7c for which A = tt, see Fig. IllB. For larger bias 7 > 7c only the d = tt 
uniform solution is stable. 

5.2.2 Bifurcation diagram 

The parameter dependence of OP solutions when interacting with OD stripes is shown in Fig. [TUA.B. 
Similar to the low order variant of this coupling energy the amplitude of the stripes pattern Ai is 
suppressed while the amplitude of the opposite mode Ai- grows. Finally both amplitudes collapse, 
leading to an orientation scotoma solution. In contrast to the low order variant this stripe pattern is stable 
for arbitrary large inter-map coupling. In case of OP rhombic solutions inter-map coupling transforms 
this solution by reducing the amplitudes A2 — A2- while increasing the amplitudes A3 ~ A^- . Without 



6A^B [4 (-4^3 + 7B^S - B5^ + S^) + Scos A {13B^ - 8BS + 6S^)] sin A = . 



(135) 




(136) 



Al 
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OD bias this solution is then transformed into the orientation scotoma stripe pattern, similar to the low 
order variant of this energy. In contrast to the low order energy, for non-zero bias the amplitudes A2 and 
A3 stay small but non-zero. 

The parameter dependence of OP solutions when interacting with OD hexagons is shown in Fig. IIOC .D. 
For a small OD bias (7 = 7*) OP rhombic solutions decay into OP stripe- like patterns. These stripe- like 
patterns stay stable also for large- inter map coupling. In case of a larger OD bias (7 = 87*), both the OP 
stripe and the OP rhombic solutions decay into the uniform PWC solution. Thus for small bias there is 
a bistability between stripe-like and uniform PWC solutions while for larger OD bias the uniform PWC 
solution is the only stable solution. The potential of OP stripe and OP rhombic solutions is shown in 
Fig. IIOE .F. In the uncoupled case as well as for small inter-map coupling strength OP stripe solutions 
are for all bias values the energetic ground state. For large inter-map coupling and a small bias (7 « 7*) 
rhombic solutions are unstable and the stripe-like solutions are energetically preferred compared to PWC 
solutions. For larger bias, however, PWC solutions are the only stable solutions for large inter-map 
coupling. 
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Figure 10. Stationary amplitudes with coupling energy U = r|z|*o^. Solid (dashed) lines: 
stable (unstable) solutions. A,B OD stripes, 7 = (blue), 7 = 7* (green), 7 = I.47* (red). C,D OD 
hexagons, 7 = 7* (blue), 7 = 87* (red). A,C Transition from OP stripe solutions, B,D Transition from 
OP rhombic solutions. E Potential, Eq. ()125p . of OP stripes and OP rhombs interacting with OD 
stripes. F Potential, Eq. (|125|) . of OP stripes, OP rhombs, and hPWC interacting with OD hexagons. 
Arrows indicate corresponding lines in the phase diagram. Fig. (jlip . 
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Figure 11. A Phase diagram with coupling energy U = ro*|z|*, Tq = 0.2, rz <C Tq. Vertical black 
lines: stability range of OD stripes, hexagons, and constant solutions. Magenta (orange) line: Stability 
border of orientation scotoma stripes. Green solid line: Stability border of rhombic solutions. Red solid 
line: Stability border of PWC solutions, red dashed line: 7c, B Course of Eq. ()136p . dashed line: A = tt. 
C Stability border between Eq. (|136p solution and the A = tt solution as a function of ro (vertical red 
line in A). 

5.2.3 Phase diagram 

The stability properties of all stationary solutions are summarized in the phase diagram Fig. 1111 Compared 
to the gradient-type interaction energy we cannot scale out the dependence on Tq. The phase diagram is 
thus plotted for ro = 0.2. We rescale the inter-map coupling strength as tB^ where B is the stationary 
amplitude of the OD hexagons. In the regime of stable OD stripes there is a transition from OP stripes 
towards the orientation scotoma stripe solution. In the regime of stable OD hexagons there is a transition 
from OP stripes towards PWC solutions (red line). The stability border of PWC solutions is strongly 
OD bias dependent and has a peak at 7 w 2j*. For small OD bias 7 the uniform solution Eq. (|136l) 
is stable. With increasing bias there is a smooth transition of this solution until at 7 = 7c the d = n 
uniform solution becomes stable. In Fig. IllC the stability border 7c between the two types of uniform 
solutions is plotted as a function of Tq. We observe that there is only a weak dependence on the control 
parameter and 7c ~ 27*. 
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5.2.4 Interaction induced pinwheel crystals 

Figure [H] illustrates the uniform solutions Eq. (|136|) for different values of the OD bias 7. For small 
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Figure 12. Bias dependent pinwheel crystals, Eq. (|136p A 7 = 7*, B 7 = 1.3 7*, C 7 = 1.6 7*, D 

7 = 27*. OP map, superimposed are the OD borders (gray), 90% ipsilateral eye dominance (black), 
and 90% contralateral eye dominance (white), = 0.2. Dashed lines mark the unit cell of the regular 
pattern. E,F Distribution of orientation preference. G Intersection angles between iso-orientation lines 
and OD borders. 



bias, the OP pattern has six pinwheels per unit cell. Two of them are located at OD maxima while one 
is located at an OD minimum. The remaining three pinwheels are located near the OD border. With 
increasing bias, these three pinwheels are pushed further away from the OD border, being attracted to 
the OD maxima. With further increasing bias three shifted pinwheels merge with the one at the OD 
maximum building a single charge 1 pinwheel centered on a contralateral peak. The remaining two 
pinwheels are located at an ispi and contra peak, respectively. Note, compared to the Braitenberg PWC 
of the A = uniform solution the charge 1 pinwheel here is located at the contralateral OD peak. Finally, 
the charge 1 pinwheels split up again into four pinwheels. With increasing bias the solution more and 
more resembles the Ipsi-center PWC (A = tt solution) which is stable also in the lower order version of 
the coupling energy. Finally, at 7/7* ~ 2 the Ipsi-center PWC becomes stable and fixed for 7 > 27*. The 
distribution of preferred orientations for different values of the bias 7 is shown in Fig.ll2E.F. reflecting the 
symmetry of each pattern. The distribution of intersection angles is shown in Fig. I12G . Remarkably, all 
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solutions show a tendency towards perpendicular intersection angles. This tendency is more pronounced 
with increasing OD bias. At about 7/7* « 1.9 parallel intersection angles are completely absent and at 
7/7* 2 there are exclusively perpendicular intersection angles. 

5.3 Gradient-type energy U = e|Vo V2;|^ 

Finally, we examine the higher order version of the gradient-type inter-map coupling. The interaction 
terms are independent of the OD shift 5. In this case the coupling strength can be rescaled as /3B^ and 
is therefore independent of the bias 7. The bias in this case only determines the stability of OD stripes, 
hexagons or the constant solution. 

5.3.1 Stationary solutions and their stability 

As for its lower order pendant a coupling to OD stripes is relatively easy to analyze. The energetic ground 
state corresponds to OP stripes with the direction perpendicular to the OD stripes for which U ~ 0. 
In addition, there are rhombic solutions with the stationary amplitudes Ai ~ Ai- = A2 = A2- = 
v/r^/(5-^80eS4). In case the OD map is a constant, o(x) = S, the gradient-type inter-map coupling 
leaves the OP unaffected. As for its lower order pendant the stationary states are therefore OP stripes 
running in an arbitrary direction and the uncoupled rhombic solutions. 

In case of OD hexagons we identified three types of non-uniform solutions. Besides stripe-like solutions of 
z(x) with one dominant mode we find rPWCs Aj = Aj- = {A, a, A) with a A and distorted rPWCs 
Aj = {Ai,A2,A3) ,Aj- = {A3,A2,Ai) with Ai ^ A2 ^ A3. Note, that distorted rPWCs arc not 
stable in case of the product-type coupling energy or the analyzed low-order coupling energies. For these 
non-uniform solutions the stationary phases are inter-map coupling strength dependent. We therefore 
calculate the stationary phases and amplitudes numerically using a Newton method and initial conditions 
close to these solutions. 

In case of OD hexagons there are further uniform solutions Aj = Aj- = A, Bj = B and ijji = "03 = 
0, 'ip2 = T^- The imaginary part of the amplitude equations, see Text SI, leads to equations for the phases 
. The ansatz Eq. (|134p leads to the stationarity condition 



(13 cos A - 5)sinA = 0. 



(138) 
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Figure 13. Stationary amplitudes with coupling energy U = e |Vz-Vo|'*, A Solid (dashed) lines: 
Stable (unstable) solutions. Blue: rPWC, green: distorted rPWC, red: hPWC. Black lines: stripe-like 
solutions. B Potential, Eq. (|125p . of OP stripes (black), OP rhombs (blue), and hPWC solutions (red). 
Arrows indicate corresponding lines in the phase diagram. Fig. (914p. 



The solutions are A = 0,A = tt, and A = ±arccos(T^) sa ±1.176 where the stationary amplitude are 



given by 



A= J (9 + ei34 (61.875 - 7.5 cos A + 4.875 cos(2A))) . 



(139) 



We calculated the stability properties of all stationary solutions by linear stability analysis considering 
perturbations of the amplitudes Aj A + a^, Aj- A + aj- and of the phases (pj — > + ipj, 



-fj - ■ This leads to a perturbation matrix AI. In general amplitude and phase perturbations 



do not decouple. We therefore calculate the eigenvalues of the perturbation M matrix numerically. It 
turns out that for this type of coupling energy only the uniform solutions with A = ±arccos(:^) are 
stable while the A = and A = tt solutions are unstable in general. 



5.3.2 Bifurcation diagram 

For increasing inter-map coupling strength the amplitudes of the OP stripe and OP rhombic solutions 
are shown in Fig. [TSlft. In case of stable OD hexagons there is a transition from rPWC (blue) towards 
distorted rPWC (green). The distorted rPWCs then decay into the hPWC (red). In case of OP stripes 
(black dashed lines) inter-map coupling leads to a slight suppression of the dominant mode and a growth 
of the remaining modes. This growth saturates at small amplitudes and thus the solution stays stripe- 
like. This stripe-like solution remains stable for arbitrary large inter-map coupling. Therefore there is a 
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bistability between hPWC solutions and stripc-likc solutions for large inter-map coupling. 
The stability borders for the rPWC and distorted rPWC solutions were obtained by calculating their 
bifurcation diagram numerically from the amplitude equations, sec Text SI. With increasing map coupling 
we observe a transition from a rPWC towards a distorted rPWC at eB"^ w 0.033 (blue dashed line in 
Fig. [M^), see also Fig. [17^. The distorted rPWC loses its stability at eB^ w 0.065 (blue sohd line in 
Fig. [Mlftl and from thereon all amplitudes are equal corresponding to the hPWC. There is a bistability 
between hPWC, rPWC, and stripe-like solutions. To calculate the inter-map coupling needed for the 
hexagonal solution to become the energetic ground state we calculated the potential Ec^. (|125l) for the 
three solutions. In case of the uniform solution Eq. (|134p the potential is given by 

V = -6A''r,-3B^ro + 27A' + ^B^ 

yt^S''e(3210- 456cosA-f 90cos(2A)) . (140) 

16 

The potential in case of the rhombic and stripe-like solutions was obtained by numerically solving the 
amplitude equations using Newtons method and initial conditions close to these solutions. Above eB^ ~ 
0.12 the hPWC is energetically preferred compared to stripe- like solutions (red dashed line in in Fig. [H]) 
and thus corresponds to the energetic ground state for large inter-map coupling, see Fig. (|13p B. 

5.3.3 Phase diagram 

We calculated the phase diagram of the coupled system in the limit ^ , shown in Fig. [TH The phase 
diagram contains the stability borders of the uncoupled OD solutions 7* , 72 , 73 , 71 ■ They correspond to 
vertical lines, as they are independent of the inter-map coupling in the limit ^ Tq. At 7 = 7* hexagons 
become stable. Stripe solutions become unstable at 7 = 73. At 7 = 73 the homogeneous solution becomes 
stable while at 7 = 74 hexagons loose their stability. In units 7/7* the borders 72 j 731 74 vary slightly 
with To , see Fig.[SJ and are drawn here for Tq = 0.2. We rescale the inter-map coupling strength as eB'^ 
where B is the stationary amplitude of the OD hexagons. The stability borders of OP solutions are then 
horizontal lines. For 7 < 7* or for 7 > 7| pinwheel free orientation stripes arc dynamically selected. For 
7* < 7 < 7I and above a critical effective coupling strength eB'^ « 0.042 hPWC solutions are stable and 
become the energetic ground state of Eq. (|125|) above eB^ w 0.117. Below eS^ w 0.065, rPWC solutions 
are stable leading to a bistability region between rPWC and hPWC solutions. We find in this region 
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Figure 14. Phase diagram with coupling energy U = e |Vz Vo|'*, for ^ Tq. Vertical lines: 
stability range of OD hexagons, green line: transition from rPWC to distorted rPWC, red line: stability 
border of hPWC, blue line: stability border of distorted rPWC. Above orange line: hPWC corresponds 
to ground state of energy. 
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Figure 15. Rhombic pinwheel crystals. A OP map with superimposed OD borders (gray), 90% 
ipsilateral eye dominance (black), and 90% contralateral eye dominance (white), 7 = 87*, Tq = 0.2. B 
Selectivity |2(x)|, white: high selectivity, black: low selectivity. 

that rhombic solutions transform into distorted rhombic solutions above an effective coupling strength of 
eS^ « 0.033. 



5.3.4 Interaction induced pinwheel crystals 

First, we studied the spatial layout of the rhombic solutions which is illustrated in Fig. [121 The rPWC 
solutions are symmetric under rotation by 180 degree. The rhombic solution has 4 pinwheels per unit cell 
and its pinwheel density is thus p = 4 cos(7r/6) ~ 3.5. One may expect that the energy term Eq. pil[) 
favors pinwheels to co-localize with OD extrema. In case of the rhombic layout there is only one pinwheel 
at an OD extremum while the other three pinwheels are located at OD saddle-points which are also 
energetically favorable positions with respect to U. The orientation selectivity \z{x)\ for the rPWC is 
shown in Fig. I15B . The pattern of selectivity is arranged in small patches of highly selective regions. 
The hexagonal layout of the two stable uniform solutions is shown in Fig. [111 The A = ± arccos(5/13) 
solutions have six pinwheels per unit cell. Their pinwheel density is therefore p = 6cos7r/6 « 5.2. Three 
pinwheels of the same topological charge are located at the extrema of the OD map. Two of these are 
located at the OD maximum while one is located at the OD minimum. The remaining three pinwheels 
are not at an OD extremum but near the OD border. The distance to the OD border depends on the 
OD bias, see Fig. I16D. For a small bias (7 sa 7*) these three pinwheels are close to the OD borders 
and with increasing bias the OD border moves away from the pinwheels. The pinwheel in the center 
of the OP hexagon is at the contralateral OD peak. Because these pinwheels organize most of the 
map while the others essentially only match one OP hexagon to its neighbors we refer to this pinwheel 
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Figure 16. Contra-center pinwheel crystals. A,B OP map, superimposed are tlie OD borders 
(gray), 90% ipsilateral eye dominance (black), and 90% contralateral eye dominance (white), 
To = 0.2,7 = 87*. A A = arccos(5/13), B A = — arccos(5/13). C Distribution of orientation preference. 
D OP map with superimposed OD map for three different values 

(7 = 7*, 7 = (74 - 7*) /2 + 7*, 7 = 74) of the OD bias. E Selectivity |2;(x)|, white: high selectivity, 
black: low selectivity. F Distribution of intersection angles. 

crystal as the Contra-center pinwheel crystal. Note, that some pinwheels are located at the vertices of 
the hexagonal pattern. Pinwheels located between these vertices (on the edge) are not in the middle 
of this edge. Solutions with A = ± arccos(5/13) are therefore not symmetric under a rotation by 60 
degree but symmetric under a rotation by 120 degree. Therefore the solution A = + arccos(5/13) cannot 
be transformed into the solution A = — arccos(5/13) by a rotation of the OD and OP pattern by 180 
degrees. This symmetry is also reflected by the distribution of preferred orientations, see Fig. I16F . Six 
orientations are slightly overrepresented. Compared to the Ipsi-center PWC, which have a 60° symmetry, 
this distribution illustrates the 120° symmetry of the pattern. The distribution of intersection angles is 
continuous, see Fig. llGC . Although there is a fixed uniform solution with varying OD bias the distribution 
of intersection angles changes. The reason for this is the bias dependent change in the OD borders, see 
Fig. I16P. For all bias values there is a tendency towards perpendicular intersection angles, although for 
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Figure 17. Inter-map coupling strength dependent pinwheel positions. OD map, 
superimposed pinwheel positions (points) for different inter-map coupling strengths, 7/7* = 3. Numbers 
label pinwheels within the unit cell (dashed lines). Blue (green, red) points: pinwheel positions for 
rPWC (distorted rPWC, hPWC) solutions. A U = e |Vz-Vo|^, using stationary amplitudes from 
Fig. ([T3| (a). Positions of distorted rPWCs move continuously (pinwheel 1,3,4). B U = r |z|^o^, using 
stationary amplitudes from Fig. (|TO| D. Positions of rPWCs move continuously (pinwheel 5,6). 

low OD bias there is an additional small peak at parallel intersection angles. The orientation selectivity 
\z{x.)\ for the hPWC is shown in Fig. [T6E. The pattern shows hexagonal bands of high selectivity. 
Finally, we study changes in pinwheel positions during the transition from a rPWC towards a hPWC i.e. 
with increasing inter-map coupling strength. In case of the higher order gradient-type coupling energy 
there is a transition towards a contra-center PWC, see Fig. [TTlft. In the regime where the distorted 
rPWC is stable, three of the four pinwheels of the rPWC are moving either from an OD saddle-point 
to a position near an OD border (pinwheel 1 and 3) or from an OD saddle-point to an OD extremum 
(pinwheel 4). One pinwheel (pinwheel 2) is fixed in space. At the transition to the hPWC two additional 
pinwheels are created, one near an OD border (pinwheel 5) and one at an OD extremum (pinwheel 6). 
We compare the inter-map coupling strength dependent pinwheel positions of the gradient-type coupling 
energy with those of the product type coupling energy, see Fig. (fT7|)B. In this case three (pinwheel 2,3,4) 
of the four rPWC pinwheels have a inter-map coupling strength independent position. The remaining 
pinwheel (pinwheel 1) with increasing inter-map coupling strength splits up into three pinwheels. While 
one of these three pinwheels (pinwheel 1) is fixed in space the remaining two pinwheels (pinwheel 5,6) 
move towards the extrema of OD. Thus for large inter-map coupling, where hPWC solutions are stable, 
all six pinwheels are located at OD extrema. 
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6 Summary 

We derived amplitude equations and analyzed ground states of the higher order inter-map coupling 
energies. We calculated local and global optima and derived corresponding phase diagrams. A main 
difference between phase diagrams for low order and high order coupling energies consists in the collapse 
of orientation selectivity above a critical coupling strength that occurs only in the low order models. In 
contrast, for the high order versions, orientation selectivity is preserved for arbitrarily strong inter-map 
coupling. In order to neglect the backrcaction on the dynamics of the modes B we assumed A <^ B <^ 1. 
Our results, however, show that for the stability of pinwheel crystals a finite amplitude B is necessary. A 
decrease in B cannot be compensated by another parameter (as it would be in case of the low order 
inter-map coupling energies). For a finite B even higher order corrections to the amplitude equations than 
those presented here can thus become significant. Such terms are we neglected in the present treatment. 
In part (II) of this study we numerically confirm our main results for the higher order inter-map coupling 
energies. 

From a practical point of view, the analyzed phase diagrams and pattern properties indicate that the 
higher order gradient-type coupling energy is the simplest and most convenient choice for constructing 
models that refiect the correlations of map layouts in the visual cortex. For this coupling, intersection 
angle statistics are reproduced well, pinwheels can be stabilized, and pattern collapse cannot occur. 
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Text S2: Supporting Information for 

Coordinated optimization of visual cortical maps 
(I) Symmetry-based analysis 

7 Amplitude equations for higher order coupling energies 

Here, we list the amplitude equations for the OP dynamics in case of the high order inter-map coupling 
energies U = e|Vz • Vo|^ and U = t|z|'*o"'. 




+(54- {B1B2B3+ B1B2B3) + Y,9^J-uv\Bu\^\B■u\^' 



u,v 




^ A,AiA,h%B,BiB,B, -Y,A^-AiA,hf.l^^B,BiBM 




A,AiAk-hfJ^_B,BiBkB,^ ^ A^-Ai-Akh^^l^_^BjBiBkB, 
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- E A,-AA,-h^^l,^_B,B,B,B,-Y^ A,-A-A,-h^;^l,_,_B,B,B,B, 



- 5h['^l^AiAjAuBu+iBu+2Bj — Sh^^^AiAjAuBu+iBu+2Bj 

j.u j.u 

- E 5hf^-^A^'Aj-AuBu+iBu+2Bj - ^ A,B„+iS„+2:b, 

- ^jAjA^- Bu+iBu+2Bj — ^i^j^u- Bu+iBu+2Bj 

j^u j.u 

- '^2Sh['^^_^_AiAj-Au-Bu+iBu+2Bj — '^^5hf^_ ^_ AiAj- A.^,- Bu+iBu+2B j 
J," J," 

u u 
u u 

- B„+2 - E 5^h^^'^ A, Au-^(u+i)Bu+2 

u u 

- 6h\^l^AjAuAu+iBu+2BjBi — Sh^^L,^Aj-AuAu+iBu+2BjBi 

j,u 

- ^ A_,-y4„- A('i(+i)-i3„+2-Bj-Bi - ^ 5hf^^ ^_ Aj- A^- A(u^iy Bu+2B jBi 

- J2 Shf^l\A,\''AuBu+iBu+2B, - J2 ^h'-^-J^r ? AuBu+iBu+2B^ 

- E \A,?K-Bu+iBu+2B, - J2 ^h^-u- \Aj- \^A^-Bu+iBu+2B, 

- E ShfA^^A,+2B'B,+^ - ^ Shf_ A'^-A,+2B]B,+^ 



where the indices are considered to be cychc i.e. j + 3 ~ j. All sums are considered to run from 1 to 
3. In the article, these amphtude equations are specified in case of OD stripes, Eq. (109), OD hexagons, 
Eq. (110), or a constant OD solution, Eq (111). 
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8 Coupling coefficients 

In the following, we list the non-zero elements of the coupling coefficients. 

Qiji ~ 9ijj — 9'iju — 24t ^^mz — Qiij — 12t 

= 48t gl, = 24r 

gijii = 6e + 12r g^jij = 33e + 48r ^ij^j = 6e + 12r = 3e + 48r 

gijju = -3e + 48r gy„„ = I.5e + I2t 

- 2 /(, = 2r 

= = = 24r = 48r 

fijii = 12t fijij = 48t + 33e fijjj = 12t + 6e = 48r + 3e 

/w"=48T-3e = I2t + 1.5e 

= + 48r /i^ = 0.756 + 24r f^j, = 0.756 + 24r 

f^l = 2l6 + 24r /i« = 9.756 + 24r h^l = 1.56 + 48r 

/*S = 10-5e+12r /i« = 4.8756 + I2r = 0.756 + 24r 

/^ll-;. = + 48r /^Ijl,^ = = 0.756 + 24r 

f^-^ - + 48r h^^^^ = = 0.756 + 24r 

^-i-^ = + 48r = = 0.756 + 24r 

/.«=24r /i|I] = 12r = 24r /i« . = I2r 

/i^^.l = 24t /i^^L = 12t /I'^.l = 6t /i^^L = 24t 

%2l tij til iji 

h^^l, = 24r = 24t h'^.'^l = 24t 

= 10.56 + I2r h^^,^ = 7.5e + 24T h^^^=2U + 24T 

C- = 8.256 + 12r = 3.756 + 12r 

'^S^-.= 7-5^ + 24^ '^*?-;-;=7.56 + 24r 
/^S-H- =15e + 48r 

/^gl,-,;- = 8.256 + 12t h'^^-j-i- = 3-756 + 12t 

= + 48r hf^l = 0.756 + 24r /ig, = 4.8756 + 12r 

= + 24r /^IJ = 24r /ig = 12r 

/lig, = 9.756 + 24r /ig, = 1.56 + 48r /i^ = 24r 

/i^^l 9.756 + 24t /I'^l ., = 1.56 + 48T 

/i^^l - - = 2l6 + 24r /I'^l = 156 + 48t 

/i^^l 0.756 + 24t h^-^l „ = 156 + 48t 
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24t 




= 156 4 
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h^'\ . 
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+ 24r 
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h 12t 
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It Jl 
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-f 24t 
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II... 
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ft . . . 


_ = iZT 


ft.. . = z4r 


1,(3) 

ft. . . 


— Z4:T 


^(3) 
ijl 
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1.(3) o.-i/l_ 

ft.. . . = oe + 4r 


h^^> . . = 
ti I J 
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a.. . 


= 21e + 


24t 


1,(3) 

ft .. . . = 

it 3 3 


= y. (be + z4r 


ft . . , — i.lDe + 4oT 


1,(3) 
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a. ._ ._ . - 
y « 3 


= 16.5e- 


h24T 


ft ^, - 


— i.06 -|- Z4T 


h . - . = 4.875e + 24t 


7 (3) 

V ] ] 


10.5e+ 12t 


y 3 I 


= 0.75eH 


h24T 


J. (3) 
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ft . , . = loe + 48r 
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ft.. ' 
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'J J 
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ft . . — z4T 
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I] I 
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fl . . 

11 1 
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, (4) 
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. = 48r 


/lO^ 

ft. . . . = ft. .' = 48T 
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ft. . . = 


I] I 


= 24t 






ft^j = 0.375e + 12t 
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- 


I, (5) 


-- 7.5e + - 


48r 




oi,(5) 
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. = . 




/^(^.) =48r 

III 11 J IJl 
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h'-'K =h'^:r =48r 






_ = 24r 




= 48t 


= I^' = 12r 






= 12t 







24t 



Note, that coupling coefficients involving the constant shift S only occur in case of the product-type 
inter-map coupling energy. 
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Text S3: Supporting Information for 

Coordinated optimization of visual cortical maps 
(I) Symmetry-based analysis 



Stability matrix 

Here, we state the stability matrices M used in the hnear stabihty analysis of the coupled amplitude 
equations. The stability matrix is defined by 



RcAi 
y IvaA,- I 



M 



Re^i- 
ImAi 



( 



Ml M2 
M3 MA 



ReAj- 
ImAi 
lu\A,- 



\ 



J 



for the uniform solutions Eq. (60). We separate the uncoupled contributions from the inter-map coupling 
contributions i.e. Ml = + A/l'^"\ The uncoupled contributions are given by 
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V 

where I denotes the 6x6 identity matrix and s = sinA,c = cosA,u = sin(2A+, 7r/6), w = sin(2A — 

7r/6), X = cos(2A + tt/Q), y ~ cos(2A ~ tt/6),w ~ sin(A — 7r/6), z = sin(A + tt/6). 

The coupling part in case of the low order product-type inter-map coupling energy is given by 
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The coupling part in case of the low order gradient-type inter-map coupling energy is given by 
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M4(^) = Ml('3),M2('3) = A/3('3) = 0. 

The stationary amplitudes A are given in Eq. (63), Eq. (68), Eq. (81), and Eq. (83). The stationary 
amplitudes B and the constant 5 are given by Eq. (94) and Eq. (116), respectively. 
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List of Figures 

1 Pinwheel annihilation, preservation, and generation in numerical simulations for 
different strengths of inter-map coupling and OD bias, Tq = 0.2, = 0.02. Color code of 
OP map with zero contours of OD map superimposed. A7 = 0, /3 = 0B7 = 0, /3 = 0.3 C 
and D 7 = 0.15,(3 = 0.3. Initial conditions identical in A - C, = 10"*, fcc,o = ^c,z liS 



89 



Stationary amplitudes with coupling energy U = a |zpo'^. Tq = 0.2. Solid (dashed) 
lines: stable (unstable) solutions to Eq. (OD stripes) and Eq. ([5T|) (OD hexagons). 
A,B OD stripes, 7 = (blue), 7 = 7* (green), 7 = I.47* (orange). C,D OD hexagons, 
7 = 1.47* (blue), 7 = 87* (red). A,C Transition from OP stripe solutions, B,D Transition 
from OP rhombic solutions. E Potential for OP stripes (red) and OP rhombs (blue) 
interacting with OD stripes, 7 = 0. F Potential for OP stripes and OP rhombs interacting 
with OD hexagons. Arrows indicate corresponding lines in the phase diagram. Fig. ([3|). . [26 
Phase diagram with the coupling energy U = ao^jzp, Tq = 0.2, <C ro- Vertical 
lines: stability range of OD stripes, hexagons, and constant solution. Magenta (orange) 
line: transition of stripes (rhombs) to the orientation scotoma solution. Blue line: stability 
border for the A = tt uniform solution (hPWC). Green line: stability line of stripe- like 
solutions. Red line: pattern solutions ceases to exist, see Eq. (|l5|) and Eq. ([57)1 . Blue 

region: stability region of hPWC, gray region: No pattern solution exists [27 

Ipsi-center pinwheel crystal. A OP map, superimposed are the OD borders (gray), 
90% ipsilateral eye dominance (black), and 90% contralateral eye dominance (white), ro = 
0.2,7 ~ 87*. Dashed lines mark the unit cell of the regular pattern. B Distribution of 

preferred orientations |28 

Stationary amplitudes with U = /3|Vz-Vop and OD hexagons. Solid (dashed) lines: 
Stable (unstable) solutions to Eq. ([55)1 . Transition from OP stripes towards the uniform 
solution (red), transition from OP rhombs towards the uniform solution (blue). Arrows 

indicate corresponding lines in the phase diagram. Fig. ([5]) [31 

Phase diagram with the coupling energy U = /3|Vz Vop, ^ Tq. Vertical black 
lines: stability range of OD stripes, hexagons, and constant solution. Blue line: stability 
border for the A = uniform solution. Green line: rhombic solutions become unstable. 
Orange line: stripe-like solutions become unstable. Red line: pattern solutions cease to 

exist, see Eq. (|5T|) . Gray region: No pattern solution exists [32 

The Braitenberg pinwheel crystal, A = uniform solution of Eq. (|184p . A OP map, 
superimposed are the OD borders (gray), 90% ipsilateral eye dominance (black), and 90% 
contralateral eye dominance (white), ro = 0.2,7 = 87*. Dashed lines mark the unit cell of 
the regular pattern. B Distribution of preferred orientations [33 
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8 OD patterns. A Stationary amplitudes of the OD dynamics. The course of Bst{l) 
Eq. (IMl) (blue), Bhexil) Eq. ^ (red), and of 5 (7) Eq. ^ (green) for = 0.2. The 
solutions are plotted in solid lines within their stability ranges. B OD map of macaque 
monkey. Adapted from |101| . C Details of B with stripe-like, patchy, and homogeneous 
layout |46 

9 The uncoupled OD dynamics. A Phase diagram of the OD model Eq. (|64|) . Dashed 
lines: stability border of hexagon solutions, solid line: stability border of stripe solution, 
gray regions: stability region of constant solution B Percentage of neurons dominated 
by the contralateral eye plotted for the three stationary solutions. Circles: numerically 
obtained values, solid lines: Cat and Chex |48 

10 Stationary amplitudes with coupling energy U = t|z|^o^. Solid (dashed) lines: 
stable (unstable) solutions. A,B OD stripes, 7 = (blue), 7 = 7* (green), 7 = I.47* (red). 
C,D OD hexagons, 7 = 7* (blue), 7 = 87* (red). A,C Transition from OP stripe solutions, 
B,D Transition from OP rhombic solutions. E Potential, Eq. (jl25p . of OP stripes and OP 
rhombs interacting with OD stripes. F Potential, Eq. (|125p . of OP stripes, OP rhombs, and 
hPWC interacting with OD hexagons. Arrows indicate corresponding lines in the phase 
diagram. Fig. ([11]) [62 

11 A Phase diagram with coupling energy U = to'*|z|^, Tq — 0.2, ^ Tq. Vertical black 
lines: stability range of OD stripes, hexagons, and constant solutions. Magenta (orange) 
line: Stability border of orientation scotoma stripes. Green solid line: Stability border of 
rhombic solutions. Red solid line: Stability border of PWC solutions, red dashed line: 7c, 
B Course of Eq. (|136p . dashed line: A = tt. C Stability border between Eq. (jl36p solution 

and the A = tt solution as a function of Tq (vertical red line in A) [63 

12 Bias dependent pinwheel crystals, Eq. ([nB| A 7 = 7*, B 7 = I.37*, C 7 = I.67*, D 
7 = 27*. OP map, superimposed are the OD borders (gray), 90% ipsilateral eye dominance 
(black), and 90% contralateral eye dominance (white), = 0.2. Dashed lines mark the 
unit cell of the regular pattern. E,F Distribution of orientation preference. G Intersection 
angles between iso-orientation lines and OD borders [64 
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13 Stationary amplitudes with coupling energy U = e |Vz Vo|^, A Solid (dashed) lines: 
Stable (unstable) solutions. Blue: rPWC, green: distorted rPWC, red: hPWC. Black lines: 
stripe-like solutions. B Potential, Eq. (|125p . of OP stripes (black), OP rhombs (blue), 
and hPWC solutions (red). Arrows indicate corresponding lines in the phase diagram. 

Fig. (^ m 

14 Phase diagram with coupling energy U = e|Vz Vo|^, for ^ Tq. Vertical lines: 
stability range of OD hexagons, green line: transition from rPWC to distorted rPWC, 
red line: stability border of hPWC, blue line: stability border of distorted rPWC. Above 
orange line: hPWC corresponds to ground state of energy [68 

15 Rhombic pinwheel crystals. A OP map with superimposed OD borders (gray), 90% 
ipsilateral eye dominance (black), and 90% contralateral eye dominance (white), 7 = 
37*,?'o — 0.2. B Selectivity |2;(x)|; white: high selectivity, black: low selectivity [68 

16 Contra-center pinwheel crystals. A,B OP map, superimposed are the OD borders 
(gray), 90% ipsilateral eye dominance (black), and 90% contralateral eye dominance (white), 
To = 0.2,7 = 37*. A A = arccos(5/13), B A = — arccos(5/13). C Distribution of ori- 
entation preference. D OP map with superimposed OD map for three different values 
(7 = 7*: 7 = (74 ~ 7*) /2 + 7*7 7 = 74) of the OD bias. E Selectivity |z(x)|, white: high 
selectivity, black: low selectivity. F Distribution of intersection angles |7C 

17 Inter-map coupling strength dependent pinwheel positions. OD map, superim- 
posed pinwheel positions (points) for different inter-map coupling strengths, 7/7* = 3. 
Numbers label pinwheels within the unit cell (dashed lines). Blue (green, red) points: pin- 
wheel positions for rPWC (distorted rPWC, hPWC) solutions. A U = e|Vz-Vo|^, using 
stationary amplitudes from Fig. (|13p (a'). Positions of distorted rPWCs move continuously 
(pinwheel 1,3,4). B U = r |z|'*o^, using stationary amplitudes from Fig. P^ D. Positions 

of rPWCs move continuously (pinwheel 5,6) |71 



